{"cells":[{"cell_type":"markdown","metadata":{"id":"79FAC74C4D88446C85DB130914CD918F","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":" # <center> Lecture8 : Posterior Inference & Prediction  </center>  \n \n ## <center> Instructor: Dr. Hu Chuan-Peng  </center>"},{"cell_type":"markdown","metadata":{"id":"9B1ED8438D2B4A1E93027F0BB775888B","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"## 使用后验分布进行统计推断  \n\n我们对贝叶斯分析中的概念、基本流程都已经有了大致的了解。我们知道通过结合先验与似然来得到后验，但更具体的后验能告诉我们什么呢？  \n\n⭐这节课我们关注怎么利用后验分布进行统计推断，包括：  \n\n* 后验估计  \n* 假设检验  \n* 后验预测  \n\n我们主要使用一个**睡眠情况调查**的例子来介绍这三部分内容。"},{"cell_type":"markdown","metadata":{"id":"960E52B2DEEC421CA09BCAF3ADC907DB","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 睡眠状况调查  \n\n假设我们想知道，通常而言，大学生11点前入睡的可能性是多少？ 我们把这一可能性记为$\\pi$，并就这一问题开展了调查。  \n\n> 我们可以试着使用贝叶斯分析来回答这个问题  \n\n1. 首先，我们需要确立先验。在开始调查之前，我们对这个入睡可能性有着某种信念，假设我们认为  \n服从  \n$$  \n\\pi \\sim Beta (4, 6)  \n$$  \n\n2. 接下来我们开始收集数据，假设对100位大学生展开调查，其中有14个人回答他们在11点前入睡.  我们的数据可以用二项分布来描述  \n\n$$  \nY|\\pi \\sim Bin (100, \\pi)  \n$$  \n\n这是我们熟悉的Beta-Binomial模型：  \n\n$$  \n\\pi \\sim Beta(\\alpha, \\beta)  \n$$  \n\n$$  \nY|\\pi \\sim Bin (n, \\pi)  \n$$  \n\n$$  \n\\pi | (Y = y) \\sim Beta (\\alpha + y, \\beta + n-y)  \n$$  \n\n3. 最后，后验可以写成：  \n\n$$  \n\\pi | (Y = 14) \\sim Beta (18, 92)  \n$$"},{"cell_type":"markdown","metadata":{"id":"EF2BD8BB662541EC90BB469396B1B123","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[]},"source":"**在贝叶斯的框架下，我们可以把后验分布看作是对$\\pi$的一种估计，即当前数据更可能在哪一种$\\pi$下出现。**  \n\n> 🤔当你看到这个后验分布时，你觉得它描述了关于$\\pi$怎样的一种信息呢？  \n> * a. 在11点前入睡的学生比例大概是16%  \n> * b. 在11点前入睡的学生比例最有可能是16%，但这个比例也可能是9%到26%中的某一个。  \n\n![imageName](https://www.bayesrulesbook.com/bookdown_files/figure-html/art-post-ch8-1.png)"},{"cell_type":"markdown","metadata":{"id":"611A3F71379342B291B6606A90D48098","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"## Posterior estimation  \n\n**选项b是更符合贝叶斯取向的回答**。 为什么？ ➡  \n\n* 我们来考虑另一种调查情况下得出的后验分布：  \n\n* 我们对$\\pi$的先验仍然为$\\pi \\sim Beta(4, 6)$，但我们只调查了10位学生，在他们之中，回答自己在11点前入睡的人数为0.  \n\n$$  \n\\pi \\sim Beta(4, 6)  \n$$  \n\n$$  \nY | \\pi  \\sim Bin(100, \\pi)  \n$$  \n\n$$  \n\\pi | (y = 0) \\sim \\text{Beta}(\\alpha + y, \\beta + n - y)  \n$$  \n\n* 这种情况下，$\\pi$的后验分布满足$\\pi | (Y = 0) \\sim Beta(4, 16)$"},{"cell_type":"markdown","metadata":{"id":"2F4AC71F465F4FF68EFFF0A6A1B17C48","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"下图展示了这两种情况下的后验分布  \n\n> 黑线表示众数，即在后验分布中最可能出现的值  \n\n![Image Name](https://cdn.kesci.com/upload/s366vtdj8a.png?imageView2/0/w/960/h/960)  \n\n* 在这两种分布下，$\\pi$最可能的取值都在16-17%左右  \n\n* 如果只关注众数，两种截然不同的情况却会导致相似的结论❌  \n\n* 但很明显，我们可以看到两图中$\\pi$的分布范围是不同的  \n\n\n### **为了同时了解$\\pi$的集中趋势和离散趋势，我们可以使用 *95%的可信区间(credible intervals)* 来描述$\\pi$**"},{"cell_type":"markdown","metadata":{"id":"B657AED2D8CC4619AB12A0E6535188EF","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 可信区间(credible intervals)  \n\n* 后验分布的2.5%百分位数和97.5%百分位数包含的区域共同组成了95%的后验可信区间  \n\n* 可信区间表示了参数出现在这个区域的概率  \n    * 比如，Beta(18,92) 95%的可信区间为(0.1,  0.24)，则我们可以说，大学生在11点前入睡的可能性落在(0.1, 0.24)的可能性是95%  \n    * 一种理解是，有95%的人认为大学生在11点前入睡的可能性很低，大概为 0.1 ~ 0.24  \n\n![Image Name](https://cdn.kesci.com/upload/s3676yx2sr.png?imageView2/0/w/960/h/960)  \n"},{"cell_type":"markdown","metadata":{"id":"83A959A7A1CF46BB802021157CCDC4F6","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 在报告可信区间时，95%是最常见的选项，但并不是唯一的选项。  \n\n* 我们也可以创建其他可信区间，比如Beta(18,92)的分布之下：  \n\t* 50% 可信区间：由25%百分位数和75%百分位数包含的区域  \n\t* 99%可信区间：由0.5%百分位数和99.5%百分位数包含的区域  \n\n![Image Name](https://cdn.kesci.com/upload/s36j8gacj.png?imageView2/0/w/960/h/960)  \n\n> 在95%CI的图中，由于后验分布是一个偏态分布，可以看到相较于前半部分，后半部分包含的值对应的概率更低。  \n> 当后验分布是偏态分布时，有时我们不以中位数为中心构建95%的可信区间，而是基于最高后验概率密度(众数)。  \n"},{"cell_type":"markdown","metadata":{"id":"8810B449DD86450694CAFD3B16585D89","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"\n### 最大后验概率密度区间(the highest density interval, HDI)  \n\n最大概率密度区间HDI与中心区间是非常接近的。但是当后验分布为非正态分布时，两者存在区别，比如：  \n\n\n![Image Name](https://cdn.kesci.com/upload/image/rjdef48vqj.png?imageView2/0/w/640/h/640)  \n\n- 此时，后验中心区间从中位数的两侧开始展开。  \n- 而最大概率密度区间HDI从y轴最高的两处开始展开。  "},{"cell_type":"markdown","metadata":{"id":"EDB88893F6224F1F87098D4D96407509","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"直觉理解：想象有一根水平的线从后验概率密度最高点往下滑，如果是95%HDI，则意味着黄色部分的面积占曲线下面积的95%  \n\n![Image Name](https://cdn.kesci.com/upload/s37y41n2yf.png?imageView2/0/w/640/h/640)  \n\n> source: https://mathematica.stackexchange.com/questions/173282/computing-credible-region-highest-posterior-density-from-empirical-distributio"},{"cell_type":"markdown","metadata":{"id":"263E9234BAB94454B07769A13A768E85","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 95% CI vs 95%HDI  \n\n![Image Name](https://cdn.kesci.com/upload/s36vflrg1f.png?imageView2/0/w/960/h/960) "},{"cell_type":"markdown","metadata":{"id":"6A0C78C8F04241649BC486BA784C001B","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"此外，非正态分布后验的平均值，中位数和众数并不相等。  \n\n- 其中众数为分布最高点对应的参数值  \n- 中位数左右两侧分布的面积各占50%  \n- 平均值的位置接近中位数，并且它容易分布形态的影响  \n\n![Image Name](https://cdn.kesci.com/upload/image/rjddhc2wo.png?imageView2/0/w/500/h/500)  \n\n> source: book--A Student’s Guide to Bayesian Statistics"},{"cell_type":"markdown","metadata":{"id":"7960A3C4240A4E38ADE0232CF93D9001","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"## Posterior hypothesis testing  \n\n>🤔假设我们在一篇报道中看到这样一种描述：在11点入睡的大学生少于20%  \n\n> 根据这个结论，我们可以对后验分布做出怎样的假设？  \n\n比如，我们可以认为，在后验分布$Beta(18,92)$中，绝大多数$\\pi$的取值都小于0.2  \n\n* 那么，我们可以考虑在后验分布$Beta(18,92)$中，$\\pi$小于0.2的概率  \n$$  \nP(\\pi < 0.2 \\; | \\; Y = 14) = \\int_{0}^{0.2}f(\\pi | y=14)d\\pi .  \n$$  \n\n(如图，可以理解为$\\pi$从0 ~ 0.2围成的曲线下面积)  \n\n![Image Name](https://cdn.kesci.com/upload/s36xa3xxpj.png?imageView2/0/w/960/h/960)  \n\n经计算可知，在$Beta(18, 92)$中，$\\pi$小于0.2的概率是89.4%"},{"cell_type":"markdown","metadata":{"id":"5488C23FDEEF41809BC0A2769E67AE0B","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 假设检验  \n\n如果我们把这个问题放到假设检验的框架下，那么这个问题可以写成：  \n\n* 虚无假设：在11点前入睡的大学生比例大于等于20%  \n\n* 备择假设：在11点前入睡的大学生比例小于20%  \n\n\n$$  \n\\begin{split}  \nH_0: & \\; \\; \\pi \\ge 0.2 \\\\  \nH_a: & \\; \\; \\pi < 0.2 \\\\  \n\\end{split}  \n$$  \n\n> 注意：这是一个单侧检验"},{"cell_type":"markdown","metadata":{"id":"52D6102B81024677AD93F3F83F2EE364","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[]},"source":"我们可以计算出这两个假设对应的后验概率：  \n\n* 备择假设：$\\pi < 0.2$  \n$$  \nP(H_a \\; | \\; Y=14) = 0.849  \n$$  \n\n* 虚无假设：$\\pi \\ge 0.2$  \n$$  \nP(H_0 \\; | \\; Y=14) = 1-P(H_a \\; | \\; Y=14) = 0.151  \n$$  \n\n将这二者相除，我们可以知道备择假设发生的可能性大约是虚无假设的6倍，这被称为**后验概率比(posterior odds)**  \n\n$$  \n\\text{posterior odds } = \\frac{P(H_a \\; | \\; Y=14)}{P(H_0 \\; | \\; Y=14)} \\approx 5.62.  \n$$  \n"},{"cell_type":"markdown","metadata":{"id":"C468B4F17AB5487B9BFA96B769D15EE4","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"我们说$\\pi$的后验概率表现了对$\\pi$信念在数据与先验之间的平均更新，那么我们可以看一下在**先验分布**中，虚无假设和备择假设发生的可能性  \n\n* 备择假设：$\\pi < 0.2$  \n$$  \nP(H_a) = \\int_0^{0.2} f(\\pi) d\\pi \\approx 0.0856  \n$$  \n\n* 虚无假设：$\\pi \\ge 0.2$  \n$$  \nP(H_0) = 1-P(H_a) \\approx 0.914  \n$$  \n\n* 先验概率比：  \n$$  \n\\text{Prior odds } = \\frac{P(H_a)}{P(H_0)} \\approx 0.093  \n$$"},{"cell_type":"markdown","metadata":{"id":"93B735546A3645C6817324F114359397","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 贝叶斯因子(Bayes Factor, BF)  \n\n> $$  \n> \\text{posterior odds } \\approx 5.62  \n> $$  \n> $$  \n> \\text{Prior odds } \\approx 0.093  \n> $$  \n\n从先验到后验，**数据**让备择假设的相对可能性发生了改变。  \n\n我们可以使用**贝叶斯因子**来量化**数据**对我们信念更新的程度：  \n\n$$  \n\\text{Bayes Factor} = \\frac{\\text{posterior odds }}{\\text{prior odds }} = \\frac{\\text{后验概率比}}{\\text{先验概率比}}  \n$$  \n\n代入我们的数据中  \n$$  \n\\text{Bayes Factor} = \\frac{5.62}{0.093} \\approx 60  \n$$  \n\n* 在这个例子中，贝叶斯因子大约是60，表明备择假设在后验分布出现的可能性大于在先验分布出现的可能性  \n\n"},{"cell_type":"markdown","metadata":{"id":"8F198ABA231B4109BF26C72D1DB7AD63","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"我们对贝叶斯因子进行一个总结：  \n\n$$  \n\\text{Bayes Factor}  \n= \\frac{\\text{posterior odds}}{\\text{prior odds}}  \n= \\frac{P(H_a | Y) / P(H_0 | Y)}{P(H_a) / P(H_0)}  \n$$  \n\n* BF=1，收集到的数据并没有改变备择假设$H_a$的相对可能性  \n* BF>1，收集到的数据增加了备择假设$H_a$的相对可能性。BF越大，表明支持备择假设$H_a$的证据越强  \n* BF<1，收集到的数据削弱了备择假设$H_a$的相对可能性。BF越小，表明支持备择假设$H_a$的证据越弱  "},{"cell_type":"markdown","metadata":{"id":"A564BB57F3494CE1A3F67F7D14CAA71C","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"|**贝叶斯因子** |**解读**|  \n|-|-|  \n|> 100|极强的证据支持*Ha*|  \n|30 ~ 100|非常强的证据支持*Ha*|  \n|10 ~ 30|较强的证据支持*Ha*|  \n|3 ~ 10|中等程度的证据支持*Ha*|  \n|1 ~ 3|较弱的证据支持*Ha*|  \n|1|没有证据|  \n|1/3 ~ 1|较弱的证据支持*H0*|  \n|1/10 ~ 1/3|中等程度的证据支持*H0*|  \n|1/30 ~ 1/10|较强的证据支持*H0*|  \n|1/100 ~ 1/30|非常强的证据支持*H0*|  \n|< 1/100|极强的证据支持*H0*|  \n\n> Source: 该表改编自胡传鹏等(2018)，源引用于Lee & Wagenmakers (2014)。"},{"cell_type":"markdown","metadata":{"id":"D4BB4530A4714AC882EBBAEF90FB9585","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[]},"source":"综上，后验概率(0.85)和贝叶斯因子(60)为备择假设提供了相当可信(fairly convincing)的证据。  \n\n* 不同于我们曾经在频率主义框架下所做的那样：根据*p*值的大小拒绝或无法拒绝零假设。  \n\n* 在贝叶斯学派中，不提倡一刀切标准。数据可以告诉我们后验概率或贝叶斯因子，从而告诉我们证据有多**强**。它只是衡量不同的假设发生的相对可能性。"},{"cell_type":"markdown","metadata":{"id":"F8632E236A4642E4AE98F63E6873852B","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**额外补充：双侧检验**  \n\n除了单侧检验，有时我们的问题需要用双侧检验来回答。  \n\n* 比如，或许我们想知道 11点前入睡的大学生比例是否为30%，此时备择假设和虚无假设可以写成：  \n$$  \n\\begin{split}  \nH_0: & \\; \\; \\pi = 0.3 \\\\  \nH_a: & \\; \\; \\pi \\ne 0.3 \\\\  \n\\end{split}  \n$$  \n\n* 但此时我们会遇到一个问题，因为$\\pi$的后验分布是连续的，我们知道对于连续型变量，单个值发生的概率为0(一条线)：  \n$$  \nP(\\pi = 0.3 | Y = 14) = \\int_{0.3}^{0.3} f(\\pi | y = 14)d\\pi = 0  \n$$  \n\n* 那么备择假设发生的概率为：  \n$$  \nP(\\pi \\ne 0.3 | Y = 14) = 1 - P(\\pi = 0.3 | Y = 14) = 1  \n$$  \n\n* 此时如果我们仍然沿用刚刚的公式计算后验概率比：  \n$$  \n\\text{Posterior odds } = \\frac{P(H_a \\; | \\; Y=14)}{P(H_0 \\; | \\; Y=14)} = \\frac{1}{0} = \\text{ nooooo!}  \n$$  \n\n* 既然在连续型变量中，考虑单个值的概率会遇到计算上的问题，我们可以引入**区间**，比如11点前入睡的大学生比例$\\pi$是否在(0.25, 0.35)这个区间内：  \n\n$$  \n\\begin{split}  \nH_0: & \\; \\; \\pi \\in (0.25, 0.35) \\\\  \nH_a: & \\; \\; \\pi \\notin (0.25, 0.35) \\\\  \n\\end{split}  \n$$  \n\n> 注：区间的选择需要结合可信区间来决定  \n\n*具体的计算步骤和单侧检验情况下一样，在这里不再展开*"},{"cell_type":"markdown","metadata":{"id":"CD4C52F6C584485097030334628232CC","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"## 后验预测 (Posterior prediction)  \n\n### 后验预测分布中的变异  \n\n在得到后验分布后，除了后验估计、假设检验之外，常见的统计推断还有后验预测。  \n\n> 后验预测：根据当前后验得出的对$\\pi$的估计，我们**预测**如果再收集20人的数据，在这新的20人中 11点之前入睡的人数是多少？  \n\n**注意**，在使用后验分布进行预测时，存在着两种变异：  \n\n1. 参数变异(**Posterior variability** in parameters)  \n\n    * 后验分布直观地体现了这一变异，它表现了不同π取值的相对可能性。  \n\n        * 比如在后验分布 $Beta(18,92)$ 中，π最可能的取值大约是0.16，但0.16并不是唯一的可能取值，我们可能会取到0.10、0.20等值  \n\n2. 采样变异(**Sampling variability** in the data)  \n\n    * 样本的抽取具有随机性"},{"cell_type":"markdown","metadata":{"id":"82896AC44B4C4E66994937D22429F383","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[]},"source":"**采样变异性的形式**  \n\n* 在新数据中，我们把可能出现的结果写为$Y' = y'$，$y'$可以为$\\{0,1,...,20\\}$中的任意一个数  \n\n* 在任意$\\pi$下，某种结果$y'$出现的可能性可以表示为:  \n\n$$  \nf(y'|\\pi) = P(Y' = y' | \\pi) = \\binom{20}{ y'} \\pi^{y'}(1-\\pi)^{20 - y'}  \n$$  \n\n* 随意选取三种不同的$\\pi$，画出每一种结果的可能性  \n\n![Image Name](https://cdn.kesci.com/upload/s36z989r7y.png?imageView2/0/w/960/h/960)  \n\n* 在每一个给定的$\\pi$下，都有不同的结果可能发生，反映了抽样过程中存在的变异性  \n    * 比如，当$\\pi=0.16$时，在当前样本中，最有可能出现的结果是有4人回答自己在11点前入睡，但由于抽样变异性，也有可能回答自己在11点前入睡的人数为0。"},{"cell_type":"markdown","metadata":{"id":"59E35F5284DA40E08F6738963F4CF3C2","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 后验预测分布  \n\n在当前例子下，我们感兴趣的是特定$y'$发生的可能性，那么我们既要考虑到给定$\\pi$下$y'$的变异性，也要考虑到$\\pi$取值的变异性:  \n$$  \nf(y'|\\pi) f(\\pi|y=14)  \n\\tag{1}  \n$$  \n\n* 遍历所有的$\\pi$，计算公式(1)的总和并进行平均，我们就可以获得$Y' = y'$发生的整体概率，得到了$Y'$的后验预测模型(**posterior predictive model** )  \n\n$$  \nf(y'|y=14) = P(Y'=y' \\; | \\; Y=y) = \\int_0^1 f(y'|\\pi) f(\\pi|y=14) d\\pi  \n$$"},{"cell_type":"markdown","metadata":{"id":"435B742B2E3F4B218DE1E0C7F88BB0A5","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[]},"source":"* 有了后验预测模型之后，我们回到最开始的问题，比如我们想知道，在这20人中，有$\\ge 5$人给出肯定回答的概率是多少？  \n\n    * 我们可以进行如下计算：  \n\n$$  \n\\begin{split}  \nP(Y' \\ge 5 | y = 14)  \n& = \\sum_{y' = 5}^{20} f(y' | y = 14) \\\\  \n& = f(y' = 5 | y = 14) + f(y' = 6 | y = 14) + \\cdots + f(y' = 20 | y = 14)\\\\  \n& = 0.233. \\\\  \n\\end{split}  \n$$"},{"cell_type":"markdown","metadata":{"id":"06D35B379E044E1F84544791F7305084","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[]},"source":"* 我们也可以计算后验预测分布的期望值：将$\\pi$与对应的后验预测值相乘并加和：  \n\n$$  \n\\begin{split}  \nE(Y' | y = 14)  \n& = \\sum_{y' = 0}^{20} y' f(y' | y = 14) \\\\  \n& = 0\\cdot f(y' = 0 | y = 14) + 1 \\cdot f(y' = 1 | y = 14) + \\cdots + 20\\cdot f(y' = 20 | y = 14)\\\\  \n& = 3.273. \\\\  \n\\end{split}  \n$$"},{"cell_type":"markdown","metadata":{"id":"BC7C3BB7220543BBB2BF3E9893D06837","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**代码练习**  \n\n后验预测分布如图所示，我们尝试使用代码来生成验预测分布  \n\n![Image Name](https://cdn.kesci.com/upload/s36z9nlnbb.png?imageView2/0/w/960/h/960)"},{"cell_type":"code","metadata":{"collapsed":false,"id":"8B36451D6A1840FD8D92AE52ED0F073A","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"import numpy as np\nimport pandas as pd\nimport scipy.stats as st\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\n\n# 定义绘图函数\ndef plot_pp(y_num, pp_sum):\n    plt.stem(y_num, pp_sum,                      #横轴为可能的结果，纵轴为结果对应的可能性\n            linefmt='black',\n            bottom=-1)\n\n    plt.xlabel(\"y'\", fontsize=15)\n    plt.xticks(np.arange(0, 21, 5))\n    plt.ylabel(\"f(y'|y=14)\",\n            fontsize=15)\n    plt.ylim(0, 0.25)\n    plt.yticks(np.arange(0, 0.21, 0.05))\n    sns.despine()","outputs":[],"execution_count":1},{"cell_type":"code","metadata":{"collapsed":false,"id":"134C0CDAFFB94E5D89CFCD66A6C49105","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"y_num = np.linspace(0, 20, 21)               # 在0-20之间生成21种可能出现的结果\ndef posterior_prediction_model(pi_i):\n  posterior = st.beta.pdf(pi_i, 18, 92)        # 计算每一个pi对应的后验分布概率值\n  likelihood = st.binom.pmf(y_num, 20, pi_i)   # 计算每一个pi下，不同y出现的可能性，一共有20个值\n  #===========================================================================\n  #          如何根据**后验分布**和**似然函数**计算后验预测模型\n  #                 请修改...中的值\n  #===========================================================================\n  pos_predict = ...\n\n  return pos_predict\n\n#===========================================================================\n#      请修改参数 pi_i 的...值\n#      输出结果的含义：在给定参数下，21个数据(对应0到20人相信学生在11点前睡)的可能性\n#===========================================================================\nposterior_prediction_model(pi_i = ...)","outputs":[],"execution_count":null},{"cell_type":"code","metadata":{"collapsed":false,"id":"476100D5B1B54A6992D0AA9DE381384C","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 计算在所有参数 pi 可能条件下不同的后验预测模型\npi = np.linspace(0,1,101)                      # 在0 - 1之间生成101个pi值\npp_list = []                                   # 创建一个空的列表，用来存储每一个pi下，不同y出现的可能性\nfor pi_i in pi:                                # 遍历pi\n    pos_predict = posterior_prediction_model(pi_i)\n    pp_list.append(pos_predict)                # 存入空列表中\n\npp_list = np.array(pp_list)                    # 转换为numpy数组\n\n# 查看后验预测模型\n# 每一行代表参数pi在0-1之间生成的101个值\n# 每一列代表21种可能出现的结果y'\npd.DataFrame(pp_list)","outputs":[],"execution_count":null},{"cell_type":"code","metadata":{"collapsed":false,"id":"8AD4F6A0DB39416A850A46924B440224","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 归一化后绘图\npp_sum = np.sum(pp_list,axis=0)             # 对所有pi下，不同y出现的可能性进行求和，axis=0表示按行对每一列进行求和\npp_sum /= np.sum(pp_sum)                    # 平均\n\nplot_pp(y_num, pp_sum)","outputs":[{"data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAkcAAAG2CAYAAAB1ZSLWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvyElEQVR4nO3de3QUdZ7//1fRJEGBZEQu6dghxBugKEq4ux3BgSg6rk6bAYddwJGLOerYIctxRXa8MKMRz+omysXhJusqiCbN6FkzQnS4tII3JOooCgiaEDoHYTQxrAZo6/eHX/pXTQfoJJ10d/J8nFPnWJ/+1KffTZ+WF5+q+pRhmqYpAAAASJI6RbsAAACAWEI4AgAAsCAcAQAAWBCOAAAALAhHAAAAFoQjAAAAC8IRAACABeGoGUzTVF1dnVgiCgCA9odw1Azff/+9UlJS9P3330e7FAAAEGGEIwAAAAvCEQAAgAXhCAAAwIJwBAAAYEE4AgAAsCAcAQAAWBCOAAAALAhHAAAAFoQjAAAAC8IRAACABeEIAADAgnAEAABgQTgCAACwIBwBAABYEI4AAAAsCEcAAAAWhCMAAAALwhEAAIAF4QgAAMCCcAQAAGBBOAIAALAgHAEAAFgQjgAAACwIRwAAABaEIwAAAAvCEQAAgAXhCAAAwIJwBAAAYEE4AgAAsCAcAQAAWBCOAAAALAhHAAAAFoQjAAAAC8IRAACABeEIAADAgnAEAABg0TnaBQCN8fv98nq98vl8stvtcjqdstls0S4LANABEI4Qczwej9xut/bv3x9oczgcKi4ulsvlimJlAICOgNNqiCkej0e5ublBwUiSqqurlZubK4/HE6XKAAAdhWGaphntIuJNXV2dUlJSVFtbq+Tk5GiX0274/X7169cvJBidYBiGHA6H9u3bxyk2AECrYeYIMcPr9Z4yGEmSaZqqqqqS1+ttw6oAAB0N4Qgxw+fzRbQfAADNQThCzLDb7RHtBwBAcxCOEDOcTqccDocMw2j0dcMwlJ6eLqfT2caVAQA6EsIRYobNZlNxcbEkhQSkE/tFRUVcjA0AaFWEI8QUl8ulkpISpaWlBbU7HA6VlJSwzhEAoNVxK38zcCt/6zvxZyxJZWVlysnJYcYIANAmmDlCTLIGoezsbIIRAKDNEI4AAAAsCEcAAAAWhCMAAAALwhEAAIAF4QgAAMCCcAQAAGBBOAIAALAgHAEAAFgQjgAAACwIRwAAABaEIwAAAAvCEQAAgAXhCAAAwIJwBAAAYEE4AgAAsCAcAQAAWBCOAAAALAhHAAAAFp2jXQDil9/vl9frlc/nk91ul9PplM1mi3ZZAAC0COEIzeLxeOR2u7V///5Am8PhUHFxsVwuVxQrAwCgZTithibzeDzKzc0NCkaSVF1drdzcXHk8nihVBgBAy8V8OFq8eLEyMzPVpUsXZWVlyev1nrKvx+PR+PHj1atXLyUnJ2vUqFFav359SL/S0lJdcsklSkpK0iWXXKJ169a15kdoV/x+v9xut0zTDHntRFt+fr78fn9blwYAQETEdDhau3at8vPzNW/ePO3YsUNOp1MTJkxQZWVlo/23bNmi8ePHq6ysTNu3b9fYsWN14403aseOHYE+27Zt06RJkzRlyhR99NFHmjJliiZOnKh33323rT5WXPN6vSEzRlamaaqqquq0IRYAgFhmmI1NAcSIESNGaMiQIVqyZEmgbeDAgbr55ptVWFgY1hiXXnqpJk2apAceeECSNGnSJNXV1emvf/1roM91112nc845R2vWrAlrzLq6OqWkpKi2tlbJyclN+ETxb82aNZo8efIZ+61evVq//e1vm/0+R44cUbdu3SRJ9fX16tq1a7PHAgCgKWJ25ujo0aPavn27cnJygtpzcnK0devWsMb46aef9P3336tHjx6Btm3btoWMee211552zIaGBtXV1QVtHZXdbo9oPwAAYk3MhqNDhw7J7/erT58+Qe19+vRRTU1NWGM88cQTOnLkiCZOnBhoq6mpafKYhYWFSklJCWzp6elN+CTti9PplMPhkGEYjb5uGIbS09PldDrbuDIAACIjZsPRCSf/JWya5in/YrZas2aNHnroIa1du1a9e/du0Zhz585VbW1tYKuqqmrCJ2hfbDabiouLJYX+OZ7YLyoqYr0jAEDcitlw1LNnT9lstpAZnYMHD4bM/Jxs7dq1mj59ul566SWNGzcu6LXU1NQmj5mUlKTk5OSgrSNzuVwqKSlRWlpaULvD4VBJSQnrHAEA4lrMhqPExERlZWWpvLw8qL28vFyjR48+5XFr1qzRbbfdptWrV+uGG24IeX3UqFEhY27YsOG0YyKUy+XSZ599FtgvKyvTvn37CEYAgLgX0ytkFxQUaMqUKRo6dKhGjRqlpUuXqrKyUnl5eZJ+Pt1VXV2t5557TtLPwWjq1KkqLi7WyJEjAzNEZ511llJSUiRJbrdb2dnZWrBggW666Sa98soreuONN/TWW29F50PGMeups+zsbE6lAQDahZidOZJ+vu2+qKhI8+fP1xVXXKEtW7aorKxMGRkZkiSfzxe05tGf//xnHT9+XHfddZfsdntgc7vdgT6jR4/Wiy++qGeffVaXX365Vq1apbVr12rEiBFt/vkAAEDsiel1jmJVR17nyKo11yJinSMAQLTE9MwRAABAWyMcAQAAWBCOAAAALAhHAAAAFoQjAAAAC8IRAACABeEIAADAgnAEAABgQTgCAACwIBwBAABYEI4AAAAsCEcAAAAWhCMAAACLztEuAGhrfr9fXq9XPp9PdrtdTqdTNpst2mUBAGIE4Qgdisfjkdvt1v79+wNtDodDxcXFcrlcUawMABArOK2GDsPj8Sg3NzcoGElSdXW1cnNz5fF4olQZACCWEI7QIfj9frndbpmmGfLaibb8/Hz5/f62Lg0AEGMIR+gQvF5vyIyRlWmaqqqqktfrbcOqAACxiHCEDsHn80W0HwCg/SIcoUOw2+0R7QcAaL8IR+gQnE6nHA6HDMNo9HXDMJSeni6n09nGlQEAYg3hCB2CzWZTcXGxJIUEpBP7RUVFrHcEACAcoeNwuVwqKSlRWlpaULvD4VBJSQnrHAEAJLEIJDoYl8ulcePGKSUlRZJUVlamnJwcZowAAAHMHKHDsQah7OxsghEAIAjhCAAAwIJwBAAAYEE4AgAAsCAcAQAAWBCOAAAALAhHAAAAFoQjAAAAC8IRAACABeEIAADAgnAEAABgQTgCAACwIBwBAABYEI4AAAAsCEcAAAAWhCMAAAALwhEAAIAF4QgAAMCCcAQAAGBBOAIAALAgHAEAAFgQjgAAACwIRwAAABaEIwAAAAvCEQAAgAXhCAAAwIJwBAAAYEE4AgAAsCAcAQAAWBCOAAAALAhHAAAAFp1bOsCxY8f0xRdf6JtvvlFtba1SUlLUq1cv9e/fXwkJCZGoEQAAoM00Kxx98803WrVqlV577TW99957amhoCOnTpUsXDR8+XDfccIOmTZumXr16tbhYAACA1maYpmmG23n37t164IEHtG7dOh09elSS1LNnT/Xv3189evRQcnKyamtr9e233+rzzz/X4cOHJUmJiYlyuVyaP3++Lrzwwtb5JG2orq5OKSkpqq2tVXJycrTLiZojR46oW7dukqT6+np17dq1w48NAIh/Yc8c/f73v9fSpUvl9/s1duxYTZ48WWPGjFFmZuYpj9m7d682btyo1atX66WXXlJpaalmzZqlp59+OiLFAwAARFrYM0dnn322Zs2apXvvvVdpaWlNfqPq6mo9/vjjWr58uY4cOdLk42MJM0c/i9fZHWaOAACnE3Y4qqmpUWpqaovfMFLjRBPh6GfxGmAIRwCA0wn7Vv5IBZp4D0YAAKB9Y50jAAAAi4iHo5kzZyorKyti4y1evFiZmZnq0qWLsrKy5PV6T9nX5/Np8uTJ6t+/vzp16qT8/PyQPqtWrZJhGCHbjz/+GLGaAQBA/Ip4ONqzZ48qKioiMtbatWuVn5+vefPmaceOHXI6nZowYYIqKysb7d/Q0KBevXpp3rx5Gjx48CnHTU5Ols/nC9q6dOkSkZoBAEB8C/tW/jvvvDOsfrt27QrpbxiGFi1a1MTSpCeffFLTp0/XjBkzJElFRUVav369lixZosLCwpD+/fr1U3FxsSRp5cqVpxzXMAyufQIAAI0KOxw988wzMgxD4dzcZhiGnnnmmaD9poajo0ePavv27brvvvuC2nNycrR169YmjXWy+vp6ZWRkyO/364orrtAf//hHXXnllafs39DQELQKeF1dXYveHwAAxK4mPT7EMAzNmDFDo0ePPmWfxx57TLt27dKzzz7bosIOHTokv9+vPn36BLX36dNHNTU1zR53wIABWrVqlS677DLV1dWpuLhYV111lT766CNddNFFjR5TWFiohx9+uNnvCQAA4kfY4WjTpk2aNWuWVqxYIcMw9Nhjj+kXv/hFSL9Vq1Zp165dmjZtWkQKNAwjaN80zZC2phg5cqRGjhwZ2L/qqqs0ZMgQPf3003rqqacaPWbu3LkqKCgI7NfV1Sk9Pb3ZNQAAgNgV9gXZ2dnZ+vjjjzV37lytWrVKAwYM0OrVq1utsJ49e8pms4XMEh08eDBkNqklOnXqpGHDhmn37t2n7JOUlKTk5OSgDQAAtE9NulstMTFRf/zjH/Xhhx/qggsu0JQpUzR+/Hjt2bMn4oUlJiYqKytL5eXlQe3l5eWnPa3XVKZpqqKiQna7PWJjAgCA+NWsW/kvueQSvf3223r66af1/vvv67LLLtP8+fN19OjRiBZXUFCg5cuXa+XKldq5c6dmz56tyspK5eXlSfr5dNfUqVODjqmoqFBFRYXq6+v1zTffqKKiQp999lng9Ycffljr16/X3r17VVFRoenTp6uioiIwJgAA6NiadEH2ye688079+te/1p133qmHHnpIq1evjmhAmjRpkg4fPqz58+fL5/Np0KBBKisrU0ZGhqSfF308ec0j611n27dv1+rVq5WRkaGvvvpKkvTdd99p1qxZqqmpUUpKiq688kpt2bJFw4cPj1jdAAAgfoX94Nkz8Xg8uueee3TgwAEZhiG/3x+JYWMSD579Wbw+HJYHzwIATqdFM0dWLpdL48eP14cffhipIQEAANpcxMKRJHXv3l1XX311JIcEAABoUxF/thoAAEA8a/Vw5HK5dMEFF7T22wAAAEREq4cjn88XuFMMAAAg1nFaDQAAwCLsC7Iff/zxZr3BgQMHmnUcAABANIQdju67775mPfC1pQ+KBQAAaEthh6PExEQdO3ZM999/vzp3Dn8FgOXLlzN7BAAA4kbYKWfw4MH64IMPlJubq8GDB4f9Bq+//jrhCAAAxI2wL8geMWKEJOmDDz5otWIAAACiLexwNHz4cJmmqXfffbdJb5Camqq+ffs2uTAAAIBoCPu0msvl0mWXXabu3bs36Q3WrVvX5KIAAACiJexwdPbZZzfpWiMAAIB4xCKQAAAAFoQjAAAAixaHoy1btmjXrl2B/V27dmnLli0tHRYAACAqWhyOxowZowULFgT2CwsLNXbs2JYOCwAAEBUROa1mmmYkhgEAAIg6rjkCAACwIBwBAABYEI4AAAAsCEcAAAAWhCMAAAALwhEAAIAF4QgAAMCCcAQAAGBBOAIAALDoHO0C0Pr8fr+8Xq98Pp/sdrucTqdsNlu0ywIAICa1OBxdffXVGjBgQGB/wIABys7ObumwiBCPxyO32639+/cH2hwOh4qLi+VyuaJYGQAAsckweTBak9XV1SklJUW1tbVKTk6Odjmn5PF4lJubG/LsO8MwJEklJSUtCkhHjhxRt27dJEn19fXq2rVr84ttJ2MDAOIf1xy1U36/X263u9GHAp9oy8/Pl9/vb+vS2i2/369NmzZpzZo12rRpE3+2ABCnmh2OqqqqIlkHIszr9QadSjuZaZqqqqqS1+ttw6raL4/Ho379+mns2LGaPHmyxo4dq379+snj8US7NABAEzU7HJ1//vn61a9+pVdeeYV/Iccgn88X0X44tROnL08Oo9XV1crNzSUgAUCcaXY46tu3r8rKyuRyuZSenq558+bpyy+/jGRtaAG73R7Rfmgcpy8BoP1pdjj68ssvVV5ert/85jf69ttvVVhYqP79+2v8+PF66aWXdOzYsUjWiSZyOp1yOByBi69PZhiG0tPT5XQ627iy9oXTlwDQ/rToguxf/vKXevHFF1VdXa0nnnhC/fv315tvvqnf/va3SktL05w5c7Rz585I1YomsNlsKi4ulqSQgHRiv6ioiPWOWojTlwDQ/kTkbrUePXpo9uzZ+vTTT/XWW29p6tSp+uGHH/Rf//VfGjRokLKzs/U///M/amhoiMTbIUwul0slJSVKS0sLanc4HC2+jR8/4/QlALQ/rbLO0e7du/XEE09o6dKl//8bGYZ69uypBx54QHfddVek37JNxcs6RyecqFeSysrKlJOTE5EZo3hdiyiSY/v9fvXr10/V1dWNXndkGIYcDof27dvHLB0AxImIrXPU0NCgF154QWPHjtWAAQO0dOlS9erVS/fee6/Wr1+v6dOnq76+Xvfcc48eeeSRSL0twmD9Szk7O5u/pCOI05cA0P60eObok08+0fLly/X888/ru+++k2maGjNmjO644w65XC4lJCQE+n799dcaOXKkEhISVFlZ2eLioyXeZo5aaxYmXmZ32mJsj8eje+65R9XV1YG29PR0FRUVcfoSAOJMs5+ttmLFCi1btkzvv/++TNNUjx495Ha7lZeXp4svvrjRYzIyMjR+/Hi98MILzS4YiEUul0vjxo1rldOXAIC21exwNHPmTEnS6NGjdccdd2jixIlKSko643GXX345q2ujXeL0JQC0D80OR3feeafy8vI0aNCgJh03Z84czZkzp7lvCwAA0KqaHY4WLlwYyToAAABiQsTuVgMAAGgPCEcAAAAWYYejTp06yWazRXSbP39+a342AACAJgv7mqPs7OxTPsS0ufr16xfR8QAAAFoq7HC0adOmViwDAAAgNnDNEQAAgAXhCAAAwIJwBAAAYBH2NUe33357SJthGFqxYkVECwIAAIimsMPRqlWrQtoIRwAAoL0JOxzt27evNesAAACICWGHo4yMjNasAwAAICY0+4LsH374IZJ1AAAAxIRmh6O0tDTdfffdqqioiGA5AAAA0dXscOT3+7V48WJlZWVp2LBhWrZsmerr6yNZGwAAQJtrdjiqqanRsmXLNGzYMG3fvl15eXmy2+2aOXOm3n333UjWCAAA0GaaHY7OPvtsTZ8+Xe+8844++eQT3X333UpKStKKFSs0evRoXX755Vq4cKG+++67CJYLAADQuiKyQvall16q4uJiHThwQM8//7yys7P16aefyu12Ky0tTVOnTpXX643EWwEAALSqiD4+JDExUZMnT9Zf/vIXud1umaapH3/8Uc8//7zGjBmjwYMH67XXXmvSmIsXL1ZmZqa6dOmirKys04Ysn8+nyZMnq3///urUqZPy8/Mb7VdaWqpLLrlESUlJuuSSS7Ru3bom1QQAANqviIajt99+W7/73e903nnnqbi4WImJiZo4caL+/Oc/65e//KX+/ve/65//+Z+1fPnysMZbu3at8vPzNW/ePO3YsUNOp1MTJkxQZWVlo/0bGhrUq1cvzZs3T4MHD260z7Zt2zRp0iRNmTJFH330kaZMmaKJEydynRQAAJAkGaZpmi0Z4PDhw3ruuee0fPlyff755zJNU+eff75mzpyp22+/Xb169Qr0fe+995STk6PevXtr165dZxx7xIgRGjJkiJYsWRJoGzhwoG6++WYVFhae9tgxY8boiiuuUFFRUVD7pEmTVFdXp7/+9a+Btuuuu07nnHOO1qxZE9ZnrqurU0pKimpra5WcnBzWMdF05MgRdevWTZJUX1+vrl27xvS48Tp2a9YMAGg7zZ45evPNN3XrrbfK4XBozpw52r17t2666Sa9/vrr2rNnj/793/89KBhJ0vDhw3XDDTeE9SiSo0ePavv27crJyQlqz8nJ0datW5tbtrZt2xYy5rXXXnvaMRsaGlRXVxe0AQCA9insx4ecbPz48ZKk9PR0zZgxQzNmzJDdbj/jcenp6XI4HGfsd+jQIfn9fvXp0yeovU+fPqqpqWle0fp5CYKmjllYWKiHH3642e8JAADiR7Nnjq6//nq9+uqr2rdvn/7whz+EFYwk6bHHHmvSQ2wNwwjaN00zpK2pmjrm3LlzVVtbG9iqqqpa9P4AACB2NXvm6H//938jWUeInj17ymazhczoHDx4MGTmpylSU1ObPGZSUpKSkpKa/Z4AACB+RPRutUhKTExUVlaWysvLg9rLy8s1evToZo87atSokDE3bNjQojEBAED7EXY4mjdvnr799tsWvdnhw4d1//33h92/oKBAy5cv18qVK7Vz507Nnj1blZWVysvLk/Tz6a6pU6cGHVNRUaGKigrV19frm2++UUVFhT777LPA6263Wxs2bNCCBQv0+eefa8GCBXrjjTdOuSYSAADoYMwwJSYmmsnJyea//du/mRUVFeEeZpqmaX7wwQem2+02u3fvbiYlJTXp2EWLFpkZGRlmYmKiOWTIEHPz5s2B16ZNm2ZeffXVQf0lhWwZGRlBfV5++WWzf//+ZkJCgjlgwACztLS0STXV1taaksza2tomHRct9fX1gT+L+vr6mB83XsduzZoBAG0n7HWO9uzZo7lz56q0tFSGYWjgwIEaM2aMhg0bpv79++ucc85R9+7dVVdXp3/84x/6/PPP9d5772njxo3as2ePTNPUb37zGxUWFur888+PdMZrU6xz1LrjxuvYrHMEAO1D2BdkX3jhhXr55Zf14YcfasmSJVq7dq0WL1582ru8TNNUt27dNGPGDN15552nXLUaAAAgVjT5brUhQ4Zo2bJlKioq0ubNm+X1evXxxx/r4MGDqq2tVUpKinr37q3BgwfL6XQqOzubf0EDAIC4EXY4uuaaa3Tdddfp3nvvlSRt375dF154oa6//vpWKw4AAKCthX232qZNm/T5558H9seMGaMFCxa0SlEAAADREnY4SkxM1JEjR4LawryWGwAAIG406YLsN998U5s3b1ZmZqakn+/IqaysDOv4vn37Nq9CAACANhR2OJo1a5by8/N1zTXXBNpKS0tVWlp6xmMNw9Dx48ebVyEAAEAbCjsc3XPPPXI4HHrllVe0f/9+bdy4Ub1799aAAQNasz4AAIA21aRb+V0ul1wulySpU6dOmjBhglauXNkqhQEAAERDk9c5OuHBBx/UlVdeGclaAAAAoq5F4QgAAKC9CftWfgAAgI6AcAQAAGBBOAIAALAgHAEAAFgQjgAAACwIRwAAABaEIwAAAAvCEQAAgAXhCAAAwIJwBAAAYEE4AgAAsCAcAQAAWBCOAAAALAhHAAAAFoQjAAAAC8IRAACABeEIAADAgnAEAABgQTgCAACwIBwBAABYEI4AAAAsCEcAAAAWhCMAAAALwhEAAIAF4QgAAMCCcAQAAGBBOAIAALAgHAEAAFh0jnYBAM7M7/fL6/XK5/PJbrfL6XTKZrNFuywAaJcIR0CM83g8crvd2r9/f6DN4XCouLhYLpcripUBQPvEaTUghnk8HuXm5gYFI0mqrq5Wbm6uPB5PlCoDgPaLcATEKL/fL7fbLdM0Q1470Zafny+/39/WpQFAu0Y4AmKU1+sNmTGyMk1TVVVV8nq9bVgVALR/hCMgRvl8voj2AwCEh3AExCi73R7RfgCA8BCOgBjldDrlcDhkGEajrxuGofT0dDmdzjauDADaN8IREKNsNpuKi4slKSQgndgvKipivSMAiDDCERDDXC6XSkpKlJaWFtTucDhUUlLCOkcA0ApYBBKIcS6XS+PGjVNKSookqaysTDk5OcwYAUArYeYIiAPWIJSdnU0wAoBWRDgCAACwIBwBAABYEI4AAAAsCEcAAAAWhCMAAAALwhEAAIAF4QgAAMCCcAQAAGBBOAIAALAgHAEAAFgQjgAAACwIRwAAABYxH44WL16szMxMdenSRVlZWfJ6vaftv3nzZmVlZalLly46//zz9cwzzwS9vmrVKhmGEbL9+OOPrfkxAABAnIjpcLR27Vrl5+dr3rx52rFjh5xOpyZMmKDKyspG++/bt0/XX3+9nE6nduzYofvvv1/33HOPSktLg/olJyfL5/MFbV26dGmLjwQAAGJc52gXcDpPPvmkpk+frhkzZkiSioqKtH79ei1ZskSFhYUh/Z955hn17dtXRUVFkqSBAwfqgw8+0H/+53/qlltuCfQzDEOpqalt8hkAAEB8idmZo6NHj2r79u3KyckJas/JydHWrVsbPWbbtm0h/a+99lp98MEHOnbsWKCtvr5eGRkZcjgc+tWvfqUdO3actpaGhgbV1dUFbQAAoH2K2XB06NAh+f1+9enTJ6i9T58+qqmpafSYmpqaRvsfP35chw4dkiQNGDBAq1at0quvvqo1a9aoS5cuuuqqq7R79+5T1lJYWKiUlJTAlp6e3sJPBwAAYlXMhqMTDMMI2jdNM6TtTP2t7SNHjtS//uu/avDgwXI6nXrppZd08cUX6+mnnz7lmHPnzlVtbW1gq6qqau7HAQAAMS5mrznq2bOnbDZbyCzRwYMHQ2aHTkhNTW20f+fOnXXuuec2ekynTp00bNiw084cJSUlKSkpqYmfAAAAxKOYnTlKTExUVlaWysvLg9rLy8s1evToRo8ZNWpUSP8NGzZo6NChSkhIaPQY0zRVUVEhu90emcIBAEBci9lwJEkFBQVavny5Vq5cqZ07d2r27NmqrKxUXl6epJ9Pd02dOjXQPy8vT19//bUKCgq0c+dOrVy5UitWrNCcOXMCfR5++GGtX79ee/fuVUVFhaZPn66KiorAmAAAoGOL2dNqkjRp0iQdPnxY8+fPl8/n06BBg1RWVqaMjAxJks/nC1rzKDMzU2VlZZo9e7YWLVqktLQ0PfXUU0G38X/33XeaNWuWampqlJKSoiuvvFJbtmzR8OHD2/zzAQCA2GOYJ65YRtjq6uqUkpKi2tpaJScnR7ucMzpy5Ii6desm6edlDLp27RrT48br2PFYMwAgVEyfVgMAAGhrhCMAAAALwhEAAIAF4QgAAMCCcAQAAGBBOAIAALCI6XWOOhK/3y+v1yufzye73S6n0ymbzRbtsgAA6HAIRzHA4/HI7XZr//79gTaHw6Hi4mK5XK4oVgYAQMfDabUo83g8ys3NDQpGklRdXa3c3Fx5PJ4oVQYAQMdEOIoiv98vt9utxhYpP9GWn58vv9/f1qUBANBhEY6iyOv1hswYWZmmqaqqKnm93jasCgCAjo1wFEU+ny+i/QAAQMsRjqLIbrdHtB8AAGg5wlEUOZ1OORwOGYbR6OuGYSg9PV1Op7ONKwMAoOMiHEWRzWZTcXGxJIUEpBP7RUVFrHcEAEAbIhxFmcvlUklJidLS0oLaHQ6HSkpKWOcIAIA2xiKQMcDlcmncuHFKSUmRJJWVlSknJ4cZIwAAooCZoxhhDULZ2dkEIwAAooSZI6AD45l+ABCKcAR0UDzTDwAax2k1oAPimX4AcGqEI6CD4Zl+AHB6hCOgg+GZfgBweoQjoIPhmX4AcHqEI6CD4Zl+AHB6hCOgg+GZfgBweoQjoIPhmX4AcHqEI6AD4pl+AHBqLAIJdFA80w8AGsfMEdCB8Uw/AAhFOAIAALAgHAEAAFgQjgAAACwIRwAAABaEIwAAAAvCEQAAgAXhCAAAwIJwBAAAYEE4AgAAsCAcAQAAWBCOAAAALAhHAAAAFoQjAAAAC8IRAACARedoFwCgffL7/fJ6vfL5fLLb7XI6nbLZbNEuCwDOiHAEIOI8Ho/cbrf2798faHM4HCouLpbL5YpiZQBwZpxWAxBRHo9Hubm5QcFIkqqrq5WbmyuPxxOlygAgPIQjABHj9/vldrtlmmbIayfa8vPz5ff727o0AAgb4QhAxHi93pAZIyvTNFVVVSWv19uGVQFA0xCOAESMz+eLaD8AiAbCEYCIsdvtEe0HANFAOAIQMU6nUw6HQ4ZhNPq6YRhKT0+X0+ls48oAIHyEIwARY7PZVFxcLEkhAenEflFREesdAYhphCMAEeVyuVRSUqK0tLSgdofDoZKSEtY5AhDzWAQSQMS5XC6NGzdOKSkpkqSysjLl5OQwYwQgLjBzBKBVWINQdnY2wQhA3GDmCEBc4ZltAFob4QhA3OCZbQDaAqfVAMQFntkGoK0QjgDEPJ7ZBqAtEY4AxLy2emab3+/Xpk2btGbNGm3atImwBXRQXHMEIOa1xTPbWvN6Ji4iB+IL4QhAzGvtZ7aduJ7p5NN2J65nasnila19EXlrBa/WDHTxODY1t83YMfMPCTPGLVq0yOzXr5+ZlJRkDhkyxNyyZctp+2/atMkcMmSImZSUZGZmZppLliwJ6VNSUmIOHDjQTExMNAcOHGh6PJ4m1VRbW2tKMmtra5t03OnU19ebkkxJZn19fcTGbc2x47Hm1hybmltv7OPHj5sOh8M0DCMwpnUzDMNMT083jx8/3uyxGxu3pWOXlpY2WrNhGKZhGGZpaWmTxzx5/JNrdzgcMTtuvI5NzfFfc1PFdDh68cUXzYSEBHPZsmXmZ599ZrrdbrNr167m119/3Wj/vXv3mmeffbbpdrvNzz77zFy2bJmZkJBglpSUBPps3brVtNls5qOPPmru3LnTfPTRR83OnTub77zzTth1EY5ad9x4HZuaW3fsE0Hj5LDR0qCxcePGUwYj67Zx48Ymjduaocs0Wy94tWagi8exqTn+a26OmA5Hw4cPN/Py8oLaBgwYYN53332N9r/33nvNAQMGBLXdcccd5siRIwP7EydONK+77rqgPtdee6156623hl0X4ah1x43Xsam59ccuLS01zzvvvKD/eaanp7fof5yrV68OKxytXr26SeO2VugyzdYLXq0Z6OJxbGqO/5qbK2avOTp69Ki2b9+u++67L6g9JydHW7dubfSYbdu2KScnJ6jt2muv1YoVK3Ts2DElJCRo27Ztmj17dkifoqKiU9bS0NCghoaGwH5tba0kqa6urikf6bSOHDkS+O+6urqI3iXTWmPHY82tOTY1t/7Y48aN0zvvvKP09HRJUklJia655hrZbLZm/x6Tk5PD7teU9/jyyy/D7jdkyJCwx5XCv3vv9ddfl9PpjPq48To2Ncd/zafSvXt3GYZx6g5tFsOaqLq62pRkvv3220HtjzzyiHnxxRc3esxFF11kPvLII0Ftb7/9tinJPHDggGmappmQkGC+8MILQX1eeOEFMzEx8ZS1PPjgg2H9C5CNjY2NjY0t9rcznfmJ2ZmjE05OdqZpnjbtNdb/5Pamjjl37lwVFBQE9n/66Sf94x//0Lnnnnv65NlEdXV1Sk9PV1VVVdj/mkVs4TuMf3yH8Y/vML61xffXvXv3074es+GoZ8+estlsqqmpCWo/ePCg+vTp0+gxqampjfbv3Lmzzj333NP2OdWYkpSUlKSkpKSgtl/84hfhfpQmS05O5gcd5/gO4x/fYfzjO4xv0fz+YnaF7MTERGVlZam8vDyovby8XKNHj270mFGjRoX037Bhg4YOHaqEhITT9jnVmAAAoGOJ2ZkjSSooKNCUKVM0dOhQjRo1SkuXLlVlZaXy8vIk/Xy6q7q6Ws8995wkKS8vTwsXLlRBQYFmzpypbdu2acWKFVqzZk1gTLfbrezsbC1YsEA33XSTXnnlFb3xxht66623ovIZAQBAbInpcDRp0iQdPnxY8+fPl8/n06BBg1RWVqaMjAxJPz8qoLKyMtA/MzNTZWVlmj17thYtWqS0tDQ99dRTuuWWWwJ9Ro8erRdffFH/8R//oT/84Q+64IILtHbtWo0YMaLNP9/JkpKS9OCDD4acwkP84DuMf3yH8Y/vML7FwvdnmGYjj7kGAADooGL2miMAAIBoIBwBAABYEI4AAAAsCEcAAAAWhKMYsnjxYmVmZqpLly7KysqS1+uNdkkI00MPPSTDMIK21NTUaJeFU9iyZYtuvPFGpaWlyTAM/eUvfwl63TRNPfTQQ0pLS9NZZ52lMWPG6NNPP41OsWjUmb7D2267LeQ3OXLkyOgUixCFhYUaNmyYunfvrt69e+vmm2/WF198EdQnmr9DwlGMWLt2rfLz8zVv3jzt2LFDTqdTEyZMCFqqALHt0ksvlc/nC2yffPJJtEvCKRw5ckSDBw/WwoULG3398ccf15NPPqmFCxfq/fffV2pqqsaPH6/vv/++jSvFqZzpO5Sk6667Lug3WVZW1oYV4nQ2b96su+66S++8847Ky8t1/Phx5eTkBD2wOqq/w9M+eQ1tZvjw4WZeXl5Q24ABA8z77rsvShWhKR588EFz8ODB0S4DzSDJXLduXWD/p59+MlNTU83HHnss0Pbjjz+aKSkp5jPPPBOFCnEmJ3+Hpmma06ZNM2+66aao1IOmO3jwoCnJ3Lx5s2ma0f8dMnMUA44ePart27crJycnqD0nJ0dbt26NUlVoqt27dystLU2ZmZm69dZbtXfv3miXhGbYt2+fampqgn6PSUlJuvrqq/k9xplNmzapd+/euvjiizVz5kwdPHgw2iXhFGprayVJPXr0kBT93yHhKAYcOnRIfr8/5OG3ffr0CXlILmLTiBEj9Nxzz2n9+vVatmyZampqNHr0aB0+fDjapaGJTvzm+D3GtwkTJuiFF17Q3/72Nz3xxBN6//33dc0116ihoSHapeEkpmmqoKBA//RP/6RBgwZJiv7vMKYfH9LRGIYRtG+aZkgbYtOECRMC/33ZZZdp1KhRuuCCC/Tf//3fKigoiGJlaC5+j/Ft0qRJgf8eNGiQhg4dqoyMDL322mtyuVxRrAwnu/vuu/Xxxx83+ozTaP0OmTmKAT179pTNZgtJwwcPHgxJzYgPXbt21WWXXabdu3dHuxQ00Ym7DPk9ti92u10ZGRn8JmPM73//e7366qvauHGjHA5HoD3av0PCUQxITExUVlaWysvLg9rLy8s1evToKFWFlmhoaNDOnTtlt9ujXQqaKDMzU6mpqUG/x6NHj2rz5s38HuPY4cOHVVVVxW8yRpimqbvvvlsej0d/+9vflJmZGfR6tH+HnFaLEQUFBZoyZYqGDh2qUaNGaenSpaqsrFReXl60S0MY5syZoxtvvFF9+/bVwYMH9ac//Ul1dXWaNm1atEtDI+rr67Vnz57A/r59+1RRUaEePXqob9++ys/P16OPPqqLLrpIF110kR599FGdffbZmjx5chSrhtXpvsMePXrooYce0i233CK73a6vvvpK999/v3r27Klf//rXUawaJ9x1111avXq1XnnlFXXv3j0wQ5SSkqKzzjpLhmFE93fY6vfDIWyLFi0yMzIyzMTERHPIkCGBWxoR+yZNmmTa7XYzISHBTEtLM10ul/npp59GuyycwsaNG01JIdu0adNM0/z5NuIHH3zQTE1NNZOSkszs7Gzzk08+iW7RCHK67/D//u//zJycHLNXr15mQkKC2bdvX3PatGlmZWVltMvG/9PYdyfJfPbZZwN9ovk7NP5fkQAAABDXHAEAAAQhHAEAAFgQjgAAACwIRwAAABaEIwAAAAvCEQAAgAXhCAAAwIJwBAAAYEE4AoAzuO2222QYhr766qtolwKgDRCOAAAALAhHAAAAFoQjAAAAC8IRgA7l/fffl2EYuuqqq07Z5+GHH5ZhGPrTn/7UhpUBiBWEIwAdyrBhw5SVlaWtW7fq008/DXn9p59+0rPPPiubzabf/e53UagQQLQRjgB0OHfccYckafny5SGvbdiwQV9//bWuv/56nXfeeZKkVatWyTRN9evXry3LBBAlhmmaZrSLAIC2dOTIEaWlpSkhIUHV1dVKSkoKvJabm6vS0lK9+uqruvHGG6NYJYBoYeYIQIfTtWtX/cu//IsOHz6sdevWBdoPHjyoV199VWlpabr++uujWCGAaCIcAeiQ8vLyJEnLli0LtK1atUrHjh3T7bffLpvNFq3SAEQZp9UAdFijRo3Su+++q927d+uCCy5Q//79tXv3bu3du5fri4AOjJkjAB1WXl6eTNPUihUrtHnzZu3atUvjx48nGAEdHDNHADqsH374Qeedd566dOmi7OxsrV27Vi+//LJyc3OjXRqAKCIcAejQ8vPzVVxcLEnq1auXqqurlZCQEOWqAEQTp9UAdGgn1jySpNtuu41gBICZIwA477zzdODAAX3xxRe6+OKLo10OgChj5ghAh7Z161YdOHBAV199NcEIgCTCEYAO7tFHH5Uk3X333VGuBECs4LQagA5n69atWrFihf7+97/rvffeU1ZWlt577z116sS/FwFInaNdAAC0tV27dmnlypXq3r27brzxRi1cuJBgBCCAmSMAAAAL/qkEAABgQTgCAACwIBwBAABYEI4AAAAsCEcAAAAWhCMAAAALwhEAAIAF4QgAAMDi/wOKXfgkeJTjHwAAAABJRU5ErkJggg==","text/plain":["<Figure size 640x480 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":10},{"cell_type":"markdown","metadata":{"id":"DACE76B5678A4EEBB738BA5ADEA13115","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"## Posterior analysis with MCMC  \n\n在之前例子中我们接触了如何通过后验分布进行后验估计、统计推断、假设检验和后验预测。  \n\n现在，我们进入实战的环节。结合之前学习 MCMC 和 PyMC 的内容，来进行对后验分布的分析。"},{"cell_type":"markdown","metadata":{"id":"4ACAC9D5289D44728DBA1374F2899338","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 定义和拟合模型  \n\n首先，定义 PyMC 模型。借用之前的例子：  \n- 先定义一个变量 $\\pi$，代表学生相信学生在11点前睡的概率。我们假设$\\pi$服从 Beta 分布，$\\alpha$ = 4, $\\beta$ = 6。  \n- 假设数据 y = 14，代表100名学生中有14名在11点前睡觉"},{"cell_type":"code","metadata":{"collapsed":false,"id":"78F6D7D7E87B469FA7E26A88A271A437","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"skip"},"tags":[],"trusted":true},"source":"# 导入 pymc 模型包，和 arviz 等分析工具 \nimport pymc as pm\nimport arviz as az\nimport seaborn as sns\nimport scipy.stats as st\nimport numpy as np\nimport matplotlib.pyplot as plt","outputs":[],"execution_count":2},{"cell_type":"code","metadata":{"collapsed":false,"id":"EADE68CDA6814B36A7689C68E12752A6","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"# STEP 1: 定义模型\nwith pm.Model() as model:\n    pi = pm.Beta('pi', alpha=4, beta=6)               # 变量 pi，代表学生相信学生在11点前睡的概率， 服从 Beta 分布\n    Y = pm.Binomial('Y', n=100, p=pi, observed=14)    # 数据 y = 14，代表100名学生中有14名在11点前睡觉","outputs":[],"execution_count":3},{"cell_type":"markdown","metadata":{"id":"A54C0C05EC164E60ABC6702418BC5CAF","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"接着，我们通过 MCMC 采样，来对后验分布进行采样。  \n- 我们设置 4 条采样链。  \n- 采样 5500 个样本，tune 500 个样本。  \n- 总共 4 * 5000 = 20000 个样本。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"AB9609D770BD4994A40B637BB108E401","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"# STEP 2: MCMC采样\nwith model:\n    trace = pm.sample(\n        chains=4,               # 设置 4 条采样链\n        draws=5000, tune=500,   # 采样 5000 个样本，tune 500 个样本\n        random_seed=84735)      # 设置随机数，以便得到相同的结果","outputs":[{"name":"stderr","output_type":"stream","text":["Auto-assigning NUTS sampler...\n","Initializing NUTS using jitter+adapt_diag...\n","Multiprocess sampling (4 chains in 4 jobs)\n","NUTS: [pi]\n"]},{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"name":"stderr","output_type":"stream","text":["Sampling 4 chains for 500 tune and 5_000 draw iterations (2_000 + 20_000 draws total) took 4 seconds.\n"]}],"execution_count":4},{"cell_type":"markdown","metadata":{"id":"680D318528FB457CA862CE8A3EB1BB19","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 模型诊断  \n\n然后，我们对模型进行可视化和诊断。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"4437693485C54BEAACFC2A50CF8A1DC5","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":true,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"az.plot_trace(trace)","outputs":[{"data":{"text/plain":["array([[<Axes: title={'center': 'pi'}>, <Axes: title={'center': 'pi'}>]],\n","      dtype=object)"]},"execution_count":5,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/483D184FF42540E8AF00B7BC58E6D9AC/s3e0q484dy.png\">"],"text/plain":["<Figure size 1200x200 with 2 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":5},{"cell_type":"markdown","metadata":{"id":"4FC75A5972384D479E2A79E5C8DC531C","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"- 图左：四条链的后验密度分布保持一致  \n- 图右：轨迹图展现了MCMC链的随机性和独立性"},{"cell_type":"code","metadata":{"collapsed":false,"id":"45105E0EEA2B40BFB35BC130C5CAB713","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"az.plot_autocorr(trace, max_lag=10, combined=True)","outputs":[{"data":{"text/plain":["<Axes: title={'center': 'pi'}>"]},"execution_count":6,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/9387C2AA006C4F2E9D5B22572263EC11/s3e0q4ulyy.png\">"],"text/plain":["<Figure size 640x480 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":6},{"cell_type":"markdown","metadata":{"id":"0748715B463E4033A5A647CD56220D78","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[]},"source":"自相关图：链表现得 \"足够 \"像一个独立样本"},{"cell_type":"markdown","metadata":{"id":"261150669E514EBE835B700AD04EBD2E","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"诊断指标显示：  \n- Rhat 接近 1，都表明模型非常稳定  \n- ess_bulk = 20,000 个采样中有 8007 个独立样本，占 40%"},{"cell_type":"code","metadata":{"collapsed":false,"id":"B6569556101F47C7A9C8A1F069D361DF","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"az.summary(trace, kind = \"diagnostics\")","outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>mcse_mean</th>\n","      <th>mcse_sd</th>\n","      <th>ess_bulk</th>\n","      <th>ess_tail</th>\n","      <th>r_hat</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>pi</th>\n","      <td>0.0</td>\n","      <td>0.0</td>\n","      <td>8007.0</td>\n","      <td>12719.0</td>\n","      <td>1.0</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["    mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat\n","pi        0.0      0.0    8007.0   12719.0    1.0"]},"execution_count":7,"metadata":{},"output_type":"execute_result"}],"execution_count":7},{"cell_type":"markdown","metadata":{"id":"3C712C3C9CDB40039E873B92C08C9D2B","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 后验估计和假设检验"},{"cell_type":"code","metadata":{"collapsed":false,"id":"56262997CC4F4B0395E7C0E2C823EE46","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"# 绘制\nsns.distplot(trace.posterior['pi'], bins=20, label=\"mcmc\")\n\nx = np.linspace(0, 0.35, 100)\ny = st.beta.pdf(x, 18, 92)\nplt.plot(x, y, label=\"true\")\n\nplt.legend()\nplt.xlabel('pi')\nplt.ylabel('Density')\nplt.title('Posterior Distribution of pi')\nplt.show()","outputs":[{"name":"stderr","output_type":"stream","text":["/tmp/ipykernel_54/2089979640.py:2: UserWarning: \n","\n","`distplot` is a deprecated function and will be removed in seaborn v0.14.0.\n","\n","Please adapt your code to use either `displot` (a figure-level function with\n","similar flexibility) or `histplot` (an axes-level function for histograms).\n","\n","For a guide to updating your code to use the new functions, please see\n","https://gist.github.com/mwaskom/de44147ed2974457ad6372750bbe5751\n","\n","  sns.distplot(trace.posterior['pi'], bins=20, label=\"mcmc\")\n"]},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/5CB886E1402D4565B0BFE988062BF9B5/s3e0q4msjo.png\">"],"text/plain":["<Figure size 640x480 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":8},{"cell_type":"markdown","metadata":{"id":"A46C60DB4F4B4B67A0F6FADACCC6C15E","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"完整的 MCMC 近似值（橙色）与实际后验模型（蓝色）非常接近。  \n\n我们可以进一步查看后验分布的均值和置信区间。  \n- hdi_2.5% 和 hdi_97.5% 报告了MCMC链值的 2.5 和 97.5 百分位数，构成了 95% 的可信区间"},{"cell_type":"code","metadata":{"collapsed":false,"id":"CE9B4FE637D840928B81DF16BF1BF222","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"# 后验分布的均值和置信区间\naz.summary(trace, kind = \"stats\", hdi_prob = 0.95)","outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>mean</th>\n","      <th>sd</th>\n","      <th>hdi_2.5%</th>\n","      <th>hdi_97.5%</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>pi</th>\n","      <td>0.164</td>\n","      <td>0.035</td>\n","      <td>0.096</td>\n","      <td>0.233</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["     mean     sd  hdi_2.5%  hdi_97.5%\n","pi  0.164  0.035     0.096      0.233"]},"execution_count":9,"metadata":{},"output_type":"execute_result"}],"execution_count":9},{"cell_type":"code","metadata":{"collapsed":false,"id":"4890BAFC9B814DD893D3050102696FB0","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"true_posterior = st.beta(18,92)\nlower_bound, upper_bound = true_posterior.interval(0.95)\nmean = true_posterior.mean()\nprint(f\"mean:\", mean)\nprint(f\"95% Lower bound:\", lower_bound)\nprint(f\"95% Upper bound:\", upper_bound)","outputs":[{"name":"stdout","output_type":"stream","text":["mean: 0.16363636363636364\n","95% Lower bound: 0.10090843809723253\n","95% Upper bound: 0.23792856204222324\n"]}],"execution_count":10},{"cell_type":"markdown","metadata":{"id":"3847683E2CA44762A95999230A4970FA","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"回到最初的问题，我们可以进一步检验：在11点前入睡的大学生比例是否大于等于20%，即  \n* 虚无假设：在11点前入睡的大学生比例大于等于20%  \n* 备择假设：在11点前入睡的大学生比例小于20%"},{"cell_type":"code","metadata":{"collapsed":false,"id":"5422D54E48114EBFAA676BB5DE0B4636","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"proportion  = (trace.posterior['pi'] > 0.2).mean().values\nprint(f\"虚无假设：在11点前入睡的大学生比例超过20%的概率: {proportion:.3f}\")\nprint(f\"备择假设：在11点前入睡的大学生比例低于20%的概率: {1-proportion:.3f}\")","outputs":[{"name":"stdout","output_type":"stream","text":["虚无假设：在11点前入睡的大学生比例超过20%的概率: 0.154\n","备择假设：在11点前入睡的大学生比例低于20%的概率: 0.846\n"]}],"execution_count":13},{"cell_type":"markdown","metadata":{"id":"40E5980F795A43529849C96A724C3A9D","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"根据**后验概率比(posterior odds)** ：  \n\n$$  \n\\text{posterior odds } = \\frac{P(H_a \\; | \\; Y=14)}{P(H_0 \\; | \\; Y=14)} \\approx 5.62.  \n$$  \n\n择假设发生的可能性大约是虚无假设的5.5倍"},{"cell_type":"code","metadata":{"collapsed":false,"id":"5EACE0329943491199D46215E36B5671","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"posterior_odds = (1-proportion)/ proportion\nprint(f\"后验概率比: {posterior_odds:.3f}\")","outputs":[{"name":"stdout","output_type":"stream","text":["后验概率比: 5.479\n"]}],"execution_count":14},{"cell_type":"markdown","metadata":{"id":"F89AF71CBA3D45CDBF23DC6D5D0345A5","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"我们可以通过贝叶斯因子来进一步验证这个假设。  \n\n$$  \n\\text{Bayes Factor} = \\frac{\\text{posterior odds }}{\\text{prior odds }} = \\frac{\\text{后验概率比}}{\\text{先验概率比}}  \n$$"},{"cell_type":"code","metadata":{"collapsed":false,"id":"19B3861293934E779AF8E173D4DFB47C","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# Define the prior distribution\nprior_alpha = 4\nprior_beta = 6\nprior_dist = st.beta(prior_alpha, prior_beta)\n\nprior_odds = prior_dist.cdf(0.2)\n\n# Compute the Bayes factor\nbayes_factor = posterior_odds/prior_odds\n\nprint(f\"Bayes Factor:{bayes_factor:.3f}\")","outputs":[{"name":"stdout","output_type":"stream","text":["Bayes Factor:63.973\n"]}],"execution_count":15},{"cell_type":"markdown","metadata":{"id":"BF09251352E74CCEA3796D1857552301","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 后验预测  \n\n最后，我们可以利用马尔可夫链值来近似计算后验预测模型 Y′  \n- 我们假设新收集的 20 位学生的数据中，有 Y’ 人在11点前睡觉。"},{"cell_type":"markdown","metadata":{"id":"126EC7C11E074C7891CBCF744B04FCE8","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[]},"source":"🤔回顾一下，后验预测模型反映了两个变异性来源：  \n  - 后验参数的变异  \n  - 参数生成数据的变异  \n\n然而我们不必如练习时自己生成后验预测，因为 pymc 提供了直接获取后验预测的方法："},{"cell_type":"code","metadata":{"collapsed":false,"id":"B7F8EF1FD78F439ABBBEA31DCA3A30AA","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"with model: \n    posterior_pred = pm.sample_posterior_predictive(trace, var_names=['Y'])\n\n# 查看结果\nposterior_pred","outputs":[{"name":"stderr","output_type":"stream","text":["Sampling: [Y]\n"]},{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"data":{"text/html":["\n","            <div>\n","              <div class='xr-header'>\n","                <div class=\"xr-obj-type\">arviz.InferenceData</div>\n","              </div>\n","              <ul class=\"xr-sections group-sections\">\n","              \n","            <li class = \"xr-section-item\">\n","                  <input id=\"idata_posterior_predictive0688e547-c761-4e1e-ac50-87d6d9572563\" class=\"xr-section-summary-in\" type=\"checkbox\">\n","                  <label for=\"idata_posterior_predictive0688e547-c761-4e1e-ac50-87d6d9572563\" class = \"xr-section-summary\">posterior_predictive</label>\n","                  <div class=\"xr-section-inline-details\"></div>\n","                  <div class=\"xr-section-details\">\n","                      <ul id=\"xr-dataset-coord-list\" class=\"xr-var-list\">\n","                          <div style=\"padding-left:2rem;\"><div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n","<defs>\n","<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n","<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n","<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n","<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n","</symbol>\n","<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n","<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n","<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n","<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n","<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n","</symbol>\n","</defs>\n","</svg>\n","<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n"," *\n"," */\n","\n",":root {\n","  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n","  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n","  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n","  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n","  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n","  --xr-background-color: var(--jp-layout-color0, white);\n","  --xr-background-color-row-even: var(--jp-layout-color1, white);\n","  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n","}\n","\n","html[theme=dark],\n","body[data-theme=dark],\n","body.vscode-dark {\n","  --xr-font-color0: rgba(255, 255, 255, 1);\n","  --xr-font-color2: rgba(255, 255, 255, 0.54);\n","  --xr-font-color3: rgba(255, 255, 255, 0.38);\n","  --xr-border-color: #1F1F1F;\n","  --xr-disabled-color: #515151;\n","  --xr-background-color: #111111;\n","  --xr-background-color-row-even: #111111;\n","  --xr-background-color-row-odd: #313131;\n","}\n","\n",".xr-wrap {\n","  display: block !important;\n","  min-width: 300px;\n","  max-width: 700px;\n","}\n","\n",".xr-text-repr-fallback {\n","  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n","  display: none;\n","}\n","\n",".xr-header {\n","  padding-top: 6px;\n","  padding-bottom: 6px;\n","  margin-bottom: 4px;\n","  border-bottom: solid 1px var(--xr-border-color);\n","}\n","\n",".xr-header > div,\n",".xr-header > ul {\n","  display: inline;\n","  margin-top: 0;\n","  margin-bottom: 0;\n","}\n","\n",".xr-obj-type,\n",".xr-array-name {\n","  margin-left: 2px;\n","  margin-right: 10px;\n","}\n","\n",".xr-obj-type {\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-sections {\n","  padding-left: 0 !important;\n","  display: grid;\n","  grid-template-columns: 150px auto auto 1fr 20px 20px;\n","}\n","\n",".xr-section-item {\n","  display: contents;\n","}\n","\n",".xr-section-item input {\n","  display: none;\n","}\n","\n",".xr-section-item input + label {\n","  color: var(--xr-disabled-color);\n","}\n","\n",".xr-section-item input:enabled + label {\n","  cursor: pointer;\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-section-item input:enabled + label:hover {\n","  color: var(--xr-font-color0);\n","}\n","\n",".xr-section-summary {\n","  grid-column: 1;\n","  color: var(--xr-font-color2);\n","  font-weight: 500;\n","}\n","\n",".xr-section-summary > span {\n","  display: inline-block;\n","  padding-left: 0.5em;\n","}\n","\n",".xr-section-summary-in:disabled + label {\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-section-summary-in + label:before {\n","  display: inline-block;\n","  content: '►';\n","  font-size: 11px;\n","  width: 15px;\n","  text-align: center;\n","}\n","\n",".xr-section-summary-in:disabled + label:before {\n","  color: var(--xr-disabled-color);\n","}\n","\n",".xr-section-summary-in:checked + label:before {\n","  content: '▼';\n","}\n","\n",".xr-section-summary-in:checked + label > span {\n","  display: none;\n","}\n","\n",".xr-section-summary,\n",".xr-section-inline-details {\n","  padding-top: 4px;\n","  padding-bottom: 4px;\n","}\n","\n",".xr-section-inline-details {\n","  grid-column: 2 / -1;\n","}\n","\n",".xr-section-details {\n","  display: none;\n","  grid-column: 1 / -1;\n","  margin-bottom: 5px;\n","}\n","\n",".xr-section-summary-in:checked ~ .xr-section-details {\n","  display: contents;\n","}\n","\n",".xr-array-wrap {\n","  grid-column: 1 / -1;\n","  display: grid;\n","  grid-template-columns: 20px auto;\n","}\n","\n",".xr-array-wrap > label {\n","  grid-column: 1;\n","  vertical-align: top;\n","}\n","\n",".xr-preview {\n","  color: var(--xr-font-color3);\n","}\n","\n",".xr-array-preview,\n",".xr-array-data {\n","  padding: 0 5px !important;\n","  grid-column: 2;\n","}\n","\n",".xr-array-data,\n",".xr-array-in:checked ~ .xr-array-preview {\n","  display: none;\n","}\n","\n",".xr-array-in:checked ~ .xr-array-data,\n",".xr-array-preview {\n","  display: inline-block;\n","}\n","\n",".xr-dim-list {\n","  display: inline-block !important;\n","  list-style: none;\n","  padding: 0 !important;\n","  margin: 0;\n","}\n","\n",".xr-dim-list li {\n","  display: inline-block;\n","  padding: 0;\n","  margin: 0;\n","}\n","\n",".xr-dim-list:before {\n","  content: '(';\n","}\n","\n",".xr-dim-list:after {\n","  content: ')';\n","}\n","\n",".xr-dim-list li:not(:last-child):after {\n","  content: ',';\n","  padding-right: 5px;\n","}\n","\n",".xr-has-index {\n","  font-weight: bold;\n","}\n","\n",".xr-var-list,\n",".xr-var-item {\n","  display: contents;\n","}\n","\n",".xr-var-item > div,\n",".xr-var-item label,\n",".xr-var-item > .xr-var-name span {\n","  background-color: var(--xr-background-color-row-even);\n","  margin-bottom: 0;\n","}\n","\n",".xr-var-item > .xr-var-name:hover span {\n","  padding-right: 5px;\n","}\n","\n",".xr-var-list > li:nth-child(odd) > div,\n",".xr-var-list > li:nth-child(odd) > label,\n",".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n","  background-color: var(--xr-background-color-row-odd);\n","}\n","\n",".xr-var-name {\n","  grid-column: 1;\n","}\n","\n",".xr-var-dims {\n","  grid-column: 2;\n","}\n","\n",".xr-var-dtype {\n","  grid-column: 3;\n","  text-align: right;\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-var-preview {\n","  grid-column: 4;\n","}\n","\n",".xr-index-preview {\n","  grid-column: 2 / 5;\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-var-name,\n",".xr-var-dims,\n",".xr-var-dtype,\n",".xr-preview,\n",".xr-attrs dt {\n","  white-space: nowrap;\n","  overflow: hidden;\n","  text-overflow: ellipsis;\n","  padding-right: 10px;\n","}\n","\n",".xr-var-name:hover,\n",".xr-var-dims:hover,\n",".xr-var-dtype:hover,\n",".xr-attrs dt:hover {\n","  overflow: visible;\n","  width: auto;\n","  z-index: 1;\n","}\n","\n",".xr-var-attrs,\n",".xr-var-data,\n",".xr-index-data {\n","  display: none;\n","  background-color: var(--xr-background-color) !important;\n","  padding-bottom: 5px !important;\n","}\n","\n",".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",".xr-var-data-in:checked ~ .xr-var-data,\n",".xr-index-data-in:checked ~ .xr-index-data {\n","  display: block;\n","}\n","\n",".xr-var-data > table {\n","  float: right;\n","}\n","\n",".xr-var-name span,\n",".xr-var-data,\n",".xr-index-name div,\n",".xr-index-data,\n",".xr-attrs {\n","  padding-left: 25px !important;\n","}\n","\n",".xr-attrs,\n",".xr-var-attrs,\n",".xr-var-data,\n",".xr-index-data {\n","  grid-column: 1 / -1;\n","}\n","\n","dl.xr-attrs {\n","  padding: 0;\n","  margin: 0;\n","  display: grid;\n","  grid-template-columns: 125px auto;\n","}\n","\n",".xr-attrs dt,\n",".xr-attrs dd {\n","  padding: 0;\n","  margin: 0;\n","  float: left;\n","  padding-right: 10px;\n","  width: auto;\n","}\n","\n",".xr-attrs dt {\n","  font-weight: normal;\n","  grid-column: 1;\n","}\n","\n",".xr-attrs dt:hover span {\n","  display: inline-block;\n","  background: var(--xr-background-color);\n","  padding-right: 10px;\n","}\n","\n",".xr-attrs dd {\n","  grid-column: 2;\n","  white-space: pre-wrap;\n","  word-break: break-all;\n","}\n","\n",".xr-icon-database,\n",".xr-icon-file-text2,\n",".xr-no-icon {\n","  display: inline-block;\n","  vertical-align: middle;\n","  width: 1em;\n","  height: 1.5em !important;\n","  stroke-width: 0;\n","  stroke: currentColor;\n","  fill: currentColor;\n","}\n","</style><pre class='xr-text-repr-fallback'>&lt;xarray.Dataset&gt;\n","Dimensions:  (chain: 4, draw: 5000)\n","Coordinates:\n","  * chain    (chain) int64 0 1 2 3\n","  * draw     (draw) int64 0 1 2 3 4 5 6 7 ... 4993 4994 4995 4996 4997 4998 4999\n","Data variables:\n","    Y        (chain, draw) int64 25 18 26 23 5 12 10 7 ... 16 17 13 5 17 14 17\n","Attributes:\n","    created_at:                 2023-10-28T07:57:46.763814\n","    arviz_version:              0.14.0\n","    inference_library:          pymc\n","    inference_library_version:  5.6.0</pre><div class='xr-wrap' style='display:none'><div class='xr-header'><div class='xr-obj-type'>xarray.Dataset</div></div><ul class='xr-sections'><li class='xr-section-item'><input id='section-96ce1ea1-1e52-4392-a59c-9f9250ff6ccb' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-96ce1ea1-1e52-4392-a59c-9f9250ff6ccb' class='xr-section-summary'  title='Expand/collapse section'>Dimensions:</label><div class='xr-section-inline-details'><ul class='xr-dim-list'><li><span class='xr-has-index'>chain</span>: 4</li><li><span class='xr-has-index'>draw</span>: 5000</li></ul></div><div class='xr-section-details'></div></li><li class='xr-section-item'><input id='section-b4a8752d-b112-41be-ab0a-c8f662bec253' class='xr-section-summary-in' type='checkbox'  checked><label for='section-b4a8752d-b112-41be-ab0a-c8f662bec253' class='xr-section-summary' >Coordinates: <span>(2)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>chain</span></div><div class='xr-var-dims'>(chain)</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>0 1 2 3</div><input id='attrs-52df315a-b6ed-443f-a847-8a1cf054f707' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-52df315a-b6ed-443f-a847-8a1cf054f707' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-64bec313-7714-4be9-b9fc-4dbbbef705eb' class='xr-var-data-in' type='checkbox'><label for='data-64bec313-7714-4be9-b9fc-4dbbbef705eb' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([0, 1, 2, 3])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>draw</span></div><div class='xr-var-dims'>(draw)</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>0 1 2 3 4 ... 4996 4997 4998 4999</div><input id='attrs-34b02b7e-b820-46cd-aab3-5a7d40071ee3' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-34b02b7e-b820-46cd-aab3-5a7d40071ee3' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-78719e3b-8076-423a-beae-a826122b0a48' class='xr-var-data-in' type='checkbox'><label for='data-78719e3b-8076-423a-beae-a826122b0a48' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([   0,    1,    2, ..., 4997, 4998, 4999])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-4351b674-8f70-4c9d-8bf7-79b2b7e29c1a' class='xr-section-summary-in' type='checkbox'  checked><label for='section-4351b674-8f70-4c9d-8bf7-79b2b7e29c1a' class='xr-section-summary' >Data variables: <span>(1)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>Y</span></div><div class='xr-var-dims'>(chain, draw)</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>25 18 26 23 5 12 ... 13 5 17 14 17</div><input id='attrs-e621424d-c1ff-463a-8933-692c2507c4b2' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-e621424d-c1ff-463a-8933-692c2507c4b2' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-fbed9e2b-5697-4602-9d06-382834a2cc6b' class='xr-var-data-in' type='checkbox'><label for='data-fbed9e2b-5697-4602-9d06-382834a2cc6b' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([[25, 18, 26, ..., 14, 14, 10],\n","       [17, 16, 16, ..., 19, 14,  9],\n","       [22, 16, 22, ..., 21, 13, 30],\n","       [18, 17, 13, ..., 17, 14, 17]])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-8f42acff-c02f-4a0a-b83a-6656c230a692' class='xr-section-summary-in' type='checkbox'  ><label for='section-8f42acff-c02f-4a0a-b83a-6656c230a692' class='xr-section-summary' >Indexes: <span>(2)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-index-name'><div>chain</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-17fff9c3-a380-4581-8984-cba24be54068' class='xr-index-data-in' type='checkbox'/><label for='index-17fff9c3-a380-4581-8984-cba24be54068' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(Int64Index([0, 1, 2, 3], dtype=&#x27;int64&#x27;, name=&#x27;chain&#x27;))</pre></div></li><li class='xr-var-item'><div class='xr-index-name'><div>draw</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-7e3185c1-eb39-422c-a8ec-383963ee46a6' class='xr-index-data-in' type='checkbox'/><label for='index-7e3185c1-eb39-422c-a8ec-383963ee46a6' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(Int64Index([   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,\n","            ...\n","            4990, 4991, 4992, 4993, 4994, 4995, 4996, 4997, 4998, 4999],\n","           dtype=&#x27;int64&#x27;, name=&#x27;draw&#x27;, length=5000))</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-43f97594-86db-4e8a-a184-1f30532aa2a8' class='xr-section-summary-in' type='checkbox'  checked><label for='section-43f97594-86db-4e8a-a184-1f30532aa2a8' class='xr-section-summary' >Attributes: <span>(4)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'><dt><span>created_at :</span></dt><dd>2023-10-28T07:57:46.763814</dd><dt><span>arviz_version :</span></dt><dd>0.14.0</dd><dt><span>inference_library :</span></dt><dd>pymc</dd><dt><span>inference_library_version :</span></dt><dd>5.6.0</dd></dl></div></li></ul></div></div><br></div>\n","                      </ul>\n","                  </div>\n","            </li>\n","            \n","            <li class = \"xr-section-item\">\n","                  <input id=\"idata_observed_datae82f8f7d-8bd7-49ac-9ba1-a85a3f17f997\" class=\"xr-section-summary-in\" type=\"checkbox\">\n","                  <label for=\"idata_observed_datae82f8f7d-8bd7-49ac-9ba1-a85a3f17f997\" class = \"xr-section-summary\">observed_data</label>\n","                  <div class=\"xr-section-inline-details\"></div>\n","                  <div class=\"xr-section-details\">\n","                      <ul id=\"xr-dataset-coord-list\" class=\"xr-var-list\">\n","                          <div style=\"padding-left:2rem;\"><div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n","<defs>\n","<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n","<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n","<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n","<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n","</symbol>\n","<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n","<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n","<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n","<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n","<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n","</symbol>\n","</defs>\n","</svg>\n","<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n"," *\n"," */\n","\n",":root {\n","  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n","  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n","  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n","  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n","  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n","  --xr-background-color: var(--jp-layout-color0, white);\n","  --xr-background-color-row-even: var(--jp-layout-color1, white);\n","  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n","}\n","\n","html[theme=dark],\n","body[data-theme=dark],\n","body.vscode-dark {\n","  --xr-font-color0: rgba(255, 255, 255, 1);\n","  --xr-font-color2: rgba(255, 255, 255, 0.54);\n","  --xr-font-color3: rgba(255, 255, 255, 0.38);\n","  --xr-border-color: #1F1F1F;\n","  --xr-disabled-color: #515151;\n","  --xr-background-color: #111111;\n","  --xr-background-color-row-even: #111111;\n","  --xr-background-color-row-odd: #313131;\n","}\n","\n",".xr-wrap {\n","  display: block !important;\n","  min-width: 300px;\n","  max-width: 700px;\n","}\n","\n",".xr-text-repr-fallback {\n","  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n","  display: none;\n","}\n","\n",".xr-header {\n","  padding-top: 6px;\n","  padding-bottom: 6px;\n","  margin-bottom: 4px;\n","  border-bottom: solid 1px var(--xr-border-color);\n","}\n","\n",".xr-header > div,\n",".xr-header > ul {\n","  display: inline;\n","  margin-top: 0;\n","  margin-bottom: 0;\n","}\n","\n",".xr-obj-type,\n",".xr-array-name {\n","  margin-left: 2px;\n","  margin-right: 10px;\n","}\n","\n",".xr-obj-type {\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-sections {\n","  padding-left: 0 !important;\n","  display: grid;\n","  grid-template-columns: 150px auto auto 1fr 20px 20px;\n","}\n","\n",".xr-section-item {\n","  display: contents;\n","}\n","\n",".xr-section-item input {\n","  display: none;\n","}\n","\n",".xr-section-item input + label {\n","  color: var(--xr-disabled-color);\n","}\n","\n",".xr-section-item input:enabled + label {\n","  cursor: pointer;\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-section-item input:enabled + label:hover {\n","  color: var(--xr-font-color0);\n","}\n","\n",".xr-section-summary {\n","  grid-column: 1;\n","  color: var(--xr-font-color2);\n","  font-weight: 500;\n","}\n","\n",".xr-section-summary > span {\n","  display: inline-block;\n","  padding-left: 0.5em;\n","}\n","\n",".xr-section-summary-in:disabled + label {\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-section-summary-in + label:before {\n","  display: inline-block;\n","  content: '►';\n","  font-size: 11px;\n","  width: 15px;\n","  text-align: center;\n","}\n","\n",".xr-section-summary-in:disabled + label:before {\n","  color: var(--xr-disabled-color);\n","}\n","\n",".xr-section-summary-in:checked + label:before {\n","  content: '▼';\n","}\n","\n",".xr-section-summary-in:checked + label > span {\n","  display: none;\n","}\n","\n",".xr-section-summary,\n",".xr-section-inline-details {\n","  padding-top: 4px;\n","  padding-bottom: 4px;\n","}\n","\n",".xr-section-inline-details {\n","  grid-column: 2 / -1;\n","}\n","\n",".xr-section-details {\n","  display: none;\n","  grid-column: 1 / -1;\n","  margin-bottom: 5px;\n","}\n","\n",".xr-section-summary-in:checked ~ .xr-section-details {\n","  display: contents;\n","}\n","\n",".xr-array-wrap {\n","  grid-column: 1 / -1;\n","  display: grid;\n","  grid-template-columns: 20px auto;\n","}\n","\n",".xr-array-wrap > label {\n","  grid-column: 1;\n","  vertical-align: top;\n","}\n","\n",".xr-preview {\n","  color: var(--xr-font-color3);\n","}\n","\n",".xr-array-preview,\n",".xr-array-data {\n","  padding: 0 5px !important;\n","  grid-column: 2;\n","}\n","\n",".xr-array-data,\n",".xr-array-in:checked ~ .xr-array-preview {\n","  display: none;\n","}\n","\n",".xr-array-in:checked ~ .xr-array-data,\n",".xr-array-preview {\n","  display: inline-block;\n","}\n","\n",".xr-dim-list {\n","  display: inline-block !important;\n","  list-style: none;\n","  padding: 0 !important;\n","  margin: 0;\n","}\n","\n",".xr-dim-list li {\n","  display: inline-block;\n","  padding: 0;\n","  margin: 0;\n","}\n","\n",".xr-dim-list:before {\n","  content: '(';\n","}\n","\n",".xr-dim-list:after {\n","  content: ')';\n","}\n","\n",".xr-dim-list li:not(:last-child):after {\n","  content: ',';\n","  padding-right: 5px;\n","}\n","\n",".xr-has-index {\n","  font-weight: bold;\n","}\n","\n",".xr-var-list,\n",".xr-var-item {\n","  display: contents;\n","}\n","\n",".xr-var-item > div,\n",".xr-var-item label,\n",".xr-var-item > .xr-var-name span {\n","  background-color: var(--xr-background-color-row-even);\n","  margin-bottom: 0;\n","}\n","\n",".xr-var-item > .xr-var-name:hover span {\n","  padding-right: 5px;\n","}\n","\n",".xr-var-list > li:nth-child(odd) > div,\n",".xr-var-list > li:nth-child(odd) > label,\n",".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n","  background-color: var(--xr-background-color-row-odd);\n","}\n","\n",".xr-var-name {\n","  grid-column: 1;\n","}\n","\n",".xr-var-dims {\n","  grid-column: 2;\n","}\n","\n",".xr-var-dtype {\n","  grid-column: 3;\n","  text-align: right;\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-var-preview {\n","  grid-column: 4;\n","}\n","\n",".xr-index-preview {\n","  grid-column: 2 / 5;\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-var-name,\n",".xr-var-dims,\n",".xr-var-dtype,\n",".xr-preview,\n",".xr-attrs dt {\n","  white-space: nowrap;\n","  overflow: hidden;\n","  text-overflow: ellipsis;\n","  padding-right: 10px;\n","}\n","\n",".xr-var-name:hover,\n",".xr-var-dims:hover,\n",".xr-var-dtype:hover,\n",".xr-attrs dt:hover {\n","  overflow: visible;\n","  width: auto;\n","  z-index: 1;\n","}\n","\n",".xr-var-attrs,\n",".xr-var-data,\n",".xr-index-data {\n","  display: none;\n","  background-color: var(--xr-background-color) !important;\n","  padding-bottom: 5px !important;\n","}\n","\n",".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",".xr-var-data-in:checked ~ .xr-var-data,\n",".xr-index-data-in:checked ~ .xr-index-data {\n","  display: block;\n","}\n","\n",".xr-var-data > table {\n","  float: right;\n","}\n","\n",".xr-var-name span,\n",".xr-var-data,\n",".xr-index-name div,\n",".xr-index-data,\n",".xr-attrs {\n","  padding-left: 25px !important;\n","}\n","\n",".xr-attrs,\n",".xr-var-attrs,\n",".xr-var-data,\n",".xr-index-data {\n","  grid-column: 1 / -1;\n","}\n","\n","dl.xr-attrs {\n","  padding: 0;\n","  margin: 0;\n","  display: grid;\n","  grid-template-columns: 125px auto;\n","}\n","\n",".xr-attrs dt,\n",".xr-attrs dd {\n","  padding: 0;\n","  margin: 0;\n","  float: left;\n","  padding-right: 10px;\n","  width: auto;\n","}\n","\n",".xr-attrs dt {\n","  font-weight: normal;\n","  grid-column: 1;\n","}\n","\n",".xr-attrs dt:hover span {\n","  display: inline-block;\n","  background: var(--xr-background-color);\n","  padding-right: 10px;\n","}\n","\n",".xr-attrs dd {\n","  grid-column: 2;\n","  white-space: pre-wrap;\n","  word-break: break-all;\n","}\n","\n",".xr-icon-database,\n",".xr-icon-file-text2,\n",".xr-no-icon {\n","  display: inline-block;\n","  vertical-align: middle;\n","  width: 1em;\n","  height: 1.5em !important;\n","  stroke-width: 0;\n","  stroke: currentColor;\n","  fill: currentColor;\n","}\n","</style><pre class='xr-text-repr-fallback'>&lt;xarray.Dataset&gt;\n","Dimensions:  (Y_dim_0: 1)\n","Coordinates:\n","  * Y_dim_0  (Y_dim_0) int64 0\n","Data variables:\n","    Y        (Y_dim_0) int64 14\n","Attributes:\n","    created_at:                 2023-10-28T07:57:46.765403\n","    arviz_version:              0.14.0\n","    inference_library:          pymc\n","    inference_library_version:  5.6.0</pre><div class='xr-wrap' style='display:none'><div class='xr-header'><div class='xr-obj-type'>xarray.Dataset</div></div><ul class='xr-sections'><li class='xr-section-item'><input id='section-4895ac6f-0f38-4d46-a304-66428983a61c' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-4895ac6f-0f38-4d46-a304-66428983a61c' class='xr-section-summary'  title='Expand/collapse section'>Dimensions:</label><div class='xr-section-inline-details'><ul class='xr-dim-list'><li><span class='xr-has-index'>Y_dim_0</span>: 1</li></ul></div><div class='xr-section-details'></div></li><li class='xr-section-item'><input id='section-26b6ead4-0089-4b5f-b854-6e3d28974398' class='xr-section-summary-in' type='checkbox'  checked><label for='section-26b6ead4-0089-4b5f-b854-6e3d28974398' class='xr-section-summary' >Coordinates: <span>(1)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>Y_dim_0</span></div><div class='xr-var-dims'>(Y_dim_0)</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>0</div><input id='attrs-0763da5f-6ee9-47d2-859e-f994bd704f9e' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-0763da5f-6ee9-47d2-859e-f994bd704f9e' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-1830e586-d299-4eac-bb97-435a6af1d029' class='xr-var-data-in' type='checkbox'><label for='data-1830e586-d299-4eac-bb97-435a6af1d029' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([0])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-8a32adda-4bce-4401-b6a3-855667be0440' class='xr-section-summary-in' type='checkbox'  checked><label for='section-8a32adda-4bce-4401-b6a3-855667be0440' class='xr-section-summary' >Data variables: <span>(1)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>Y</span></div><div class='xr-var-dims'>(Y_dim_0)</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>14</div><input id='attrs-5fd70e38-6ab7-4bf9-9341-5e196ebfe094' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-5fd70e38-6ab7-4bf9-9341-5e196ebfe094' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-867f41ff-7a44-4a00-aacd-cf0295b0f751' class='xr-var-data-in' type='checkbox'><label for='data-867f41ff-7a44-4a00-aacd-cf0295b0f751' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([14])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-e6676c91-63d3-4f0f-aa1f-318e9906e29d' class='xr-section-summary-in' type='checkbox'  ><label for='section-e6676c91-63d3-4f0f-aa1f-318e9906e29d' class='xr-section-summary' >Indexes: <span>(1)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-index-name'><div>Y_dim_0</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-13b37e63-0c3c-4281-97c0-c869e88a37fd' class='xr-index-data-in' type='checkbox'/><label for='index-13b37e63-0c3c-4281-97c0-c869e88a37fd' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(Int64Index([0], dtype=&#x27;int64&#x27;, name=&#x27;Y_dim_0&#x27;))</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-95c5cb90-f53c-4190-96bf-8238aae07a86' class='xr-section-summary-in' type='checkbox'  checked><label for='section-95c5cb90-f53c-4190-96bf-8238aae07a86' class='xr-section-summary' >Attributes: <span>(4)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'><dt><span>created_at :</span></dt><dd>2023-10-28T07:57:46.765403</dd><dt><span>arviz_version :</span></dt><dd>0.14.0</dd><dt><span>inference_library :</span></dt><dd>pymc</dd><dt><span>inference_library_version :</span></dt><dd>5.6.0</dd></dl></div></li></ul></div></div><br></div>\n","                      </ul>\n","                  </div>\n","            </li>\n","            \n","              </ul>\n","            </div>\n","            <style> /* CSS stylesheet for displaying InferenceData objects in jupyterlab.\n"," *\n"," */\n","\n",":root {\n","  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n","  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n","  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n","  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n","  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n","  --xr-background-color: var(--jp-layout-color0, white);\n","  --xr-background-color-row-even: var(--jp-layout-color1, white);\n","  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n","}\n","\n","html[theme=dark],\n","body.vscode-dark {\n","  --xr-font-color0: rgba(255, 255, 255, 1);\n","  --xr-font-color2: rgba(255, 255, 255, 0.54);\n","  --xr-font-color3: rgba(255, 255, 255, 0.38);\n","  --xr-border-color: #1F1F1F;\n","  --xr-disabled-color: #515151;\n","  --xr-background-color: #111111;\n","  --xr-background-color-row-even: #111111;\n","  --xr-background-color-row-odd: #313131;\n","}\n","\n",".xr-wrap {\n","  display: block;\n","  min-width: 300px;\n","  max-width: 700px;\n","}\n","\n",".xr-text-repr-fallback {\n","  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n","  display: none;\n","}\n","\n",".xr-header {\n","  padding-top: 6px;\n","  padding-bottom: 6px;\n","  margin-bottom: 4px;\n","  border-bottom: solid 1px var(--xr-border-color);\n","}\n","\n",".xr-header > div,\n",".xr-header > ul {\n","  display: inline;\n","  margin-top: 0;\n","  margin-bottom: 0;\n","}\n","\n",".xr-obj-type,\n",".xr-array-name {\n","  margin-left: 2px;\n","  margin-right: 10px;\n","}\n","\n",".xr-obj-type {\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-sections {\n","  padding-left: 0 !important;\n","  display: grid;\n","  grid-template-columns: 150px auto auto 1fr 20px 20px;\n","}\n","\n",".xr-sections.group-sections {\n","  grid-template-columns: auto;\n","}\n","\n",".xr-section-item {\n","  display: contents;\n","}\n","\n",".xr-section-item input {\n","  display: none;\n","}\n","\n",".xr-section-item input + label {\n","  color: var(--xr-disabled-color);\n","}\n","\n",".xr-section-item input:enabled + label {\n","  cursor: pointer;\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-section-item input:enabled + label:hover {\n","  color: var(--xr-font-color0);\n","}\n","\n",".xr-section-summary {\n","  grid-column: 1;\n","  color: var(--xr-font-color2);\n","  font-weight: 500;\n","}\n","\n",".xr-section-summary > span {\n","  display: inline-block;\n","  padding-left: 0.5em;\n","}\n","\n",".xr-section-summary-in:disabled + label {\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-section-summary-in + label:before {\n","  display: inline-block;\n","  content: '►';\n","  font-size: 11px;\n","  width: 15px;\n","  text-align: center;\n","}\n","\n",".xr-section-summary-in:disabled + label:before {\n","  color: var(--xr-disabled-color);\n","}\n","\n",".xr-section-summary-in:checked + label:before {\n","  content: '▼';\n","}\n","\n",".xr-section-summary-in:checked + label > span {\n","  display: none;\n","}\n","\n",".xr-section-summary,\n",".xr-section-inline-details {\n","  padding-top: 4px;\n","  padding-bottom: 4px;\n","}\n","\n",".xr-section-inline-details {\n","  grid-column: 2 / -1;\n","}\n","\n",".xr-section-details {\n","  display: none;\n","  grid-column: 1 / -1;\n","  margin-bottom: 5px;\n","}\n","\n",".xr-section-summary-in:checked ~ .xr-section-details {\n","  display: contents;\n","}\n","\n",".xr-array-wrap {\n","  grid-column: 1 / -1;\n","  display: grid;\n","  grid-template-columns: 20px auto;\n","}\n","\n",".xr-array-wrap > label {\n","  grid-column: 1;\n","  vertical-align: top;\n","}\n","\n",".xr-preview {\n","  color: var(--xr-font-color3);\n","}\n","\n",".xr-array-preview,\n",".xr-array-data {\n","  padding: 0 5px !important;\n","  grid-column: 2;\n","}\n","\n",".xr-array-data,\n",".xr-array-in:checked ~ .xr-array-preview {\n","  display: none;\n","}\n","\n",".xr-array-in:checked ~ .xr-array-data,\n",".xr-array-preview {\n","  display: inline-block;\n","}\n","\n",".xr-dim-list {\n","  display: inline-block !important;\n","  list-style: none;\n","  padding: 0 !important;\n","  margin: 0;\n","}\n","\n",".xr-dim-list li {\n","  display: inline-block;\n","  padding: 0;\n","  margin: 0;\n","}\n","\n",".xr-dim-list:before {\n","  content: '(';\n","}\n","\n",".xr-dim-list:after {\n","  content: ')';\n","}\n","\n",".xr-dim-list li:not(:last-child):after {\n","  content: ',';\n","  padding-right: 5px;\n","}\n","\n",".xr-has-index {\n","  font-weight: bold;\n","}\n","\n",".xr-var-list,\n",".xr-var-item {\n","  display: contents;\n","}\n","\n",".xr-var-item > div,\n",".xr-var-item label,\n",".xr-var-item > .xr-var-name span {\n","  background-color: var(--xr-background-color-row-even);\n","  margin-bottom: 0;\n","}\n","\n",".xr-var-item > .xr-var-name:hover span {\n","  padding-right: 5px;\n","}\n","\n",".xr-var-list > li:nth-child(odd) > div,\n",".xr-var-list > li:nth-child(odd) > label,\n",".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n","  background-color: var(--xr-background-color-row-odd);\n","}\n","\n",".xr-var-name {\n","  grid-column: 1;\n","}\n","\n",".xr-var-dims {\n","  grid-column: 2;\n","}\n","\n",".xr-var-dtype {\n","  grid-column: 3;\n","  text-align: right;\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-var-preview {\n","  grid-column: 4;\n","}\n","\n",".xr-var-name,\n",".xr-var-dims,\n",".xr-var-dtype,\n",".xr-preview,\n",".xr-attrs dt {\n","  white-space: nowrap;\n","  overflow: hidden;\n","  text-overflow: ellipsis;\n","  padding-right: 10px;\n","}\n","\n",".xr-var-name:hover,\n",".xr-var-dims:hover,\n",".xr-var-dtype:hover,\n",".xr-attrs dt:hover {\n","  overflow: visible;\n","  width: auto;\n","  z-index: 1;\n","}\n","\n",".xr-var-attrs,\n",".xr-var-data {\n","  display: none;\n","  background-color: var(--xr-background-color) !important;\n","  padding-bottom: 5px !important;\n","}\n","\n",".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",".xr-var-data-in:checked ~ .xr-var-data {\n","  display: block;\n","}\n","\n",".xr-var-data > table {\n","  float: right;\n","}\n","\n",".xr-var-name span,\n",".xr-var-data,\n",".xr-attrs {\n","  padding-left: 25px !important;\n","}\n","\n",".xr-attrs,\n",".xr-var-attrs,\n",".xr-var-data {\n","  grid-column: 1 / -1;\n","}\n","\n","dl.xr-attrs {\n","  padding: 0;\n","  margin: 0;\n","  display: grid;\n","  grid-template-columns: 125px auto;\n","}\n","\n",".xr-attrs dt, dd {\n","  padding: 0;\n","  margin: 0;\n","  float: left;\n","  padding-right: 10px;\n","  width: auto;\n","}\n","\n",".xr-attrs dt {\n","  font-weight: normal;\n","  grid-column: 1;\n","}\n","\n",".xr-attrs dt:hover span {\n","  display: inline-block;\n","  background: var(--xr-background-color);\n","  padding-right: 10px;\n","}\n","\n",".xr-attrs dd {\n","  grid-column: 2;\n","  white-space: pre-wrap;\n","  word-break: break-all;\n","}\n","\n",".xr-icon-database,\n",".xr-icon-file-text2 {\n","  display: inline-block;\n","  vertical-align: middle;\n","  width: 1em;\n","  height: 1.5em !important;\n","  stroke-width: 0;\n","  stroke: currentColor;\n","  fill: currentColor;\n","}\n",".xr-wrap{width:700px!important;} </style>"],"text/plain":["Inference data with groups:\n","\t> posterior_predictive\n","\t> observed_data"]},"execution_count":16,"metadata":{},"output_type":"execute_result"}],"execution_count":16},{"cell_type":"markdown","metadata":{"id":"E7626AF71A664675B47D530AF508E7CE","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[]},"source":"返回的 `posterior_pred` 包含之前的观测数据 `observed_data`，以及预测的数据 `posterior_predictive`。  \n\n我们可以进一步观察预测结果，并进行绘图。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"533C7566A8B543CBA09D03C54270FCE1","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"observed_data = posterior_pred.observed_data['Y']\nposterior_pred = posterior_pred.posterior_predictive.to_dataframe()\nposterior_pred","outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th></th>\n","      <th>Y</th>\n","    </tr>\n","    <tr>\n","      <th>chain</th>\n","      <th>draw</th>\n","      <th></th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th rowspan=\"5\" valign=\"top\">0</th>\n","      <th>0</th>\n","      <td>25</td>\n","    </tr>\n","    <tr>\n","      <th>1</th>\n","      <td>18</td>\n","    </tr>\n","    <tr>\n","      <th>2</th>\n","      <td>26</td>\n","    </tr>\n","    <tr>\n","      <th>3</th>\n","      <td>23</td>\n","    </tr>\n","    <tr>\n","      <th>4</th>\n","      <td>5</td>\n","    </tr>\n","    <tr>\n","      <th>...</th>\n","      <th>...</th>\n","      <td>...</td>\n","    </tr>\n","    <tr>\n","      <th rowspan=\"5\" valign=\"top\">3</th>\n","      <th>4995</th>\n","      <td>13</td>\n","    </tr>\n","    <tr>\n","      <th>4996</th>\n","      <td>5</td>\n","    </tr>\n","    <tr>\n","      <th>4997</th>\n","      <td>17</td>\n","    </tr>\n","    <tr>\n","      <th>4998</th>\n","      <td>14</td>\n","    </tr>\n","    <tr>\n","      <th>4999</th>\n","      <td>17</td>\n","    </tr>\n","  </tbody>\n","</table>\n","<p>20000 rows × 1 columns</p>\n","</div>"],"text/plain":["             Y\n","chain draw    \n","0     0     25\n","      1     18\n","      2     26\n","      3     23\n","      4      5\n","...         ..\n","3     4995  13\n","      4996   5\n","      4997  17\n","      4998  14\n","      4999  17\n","\n","[20000 rows x 1 columns]"]},"execution_count":17,"metadata":{},"output_type":"execute_result"}],"execution_count":17},{"cell_type":"markdown","metadata":{"id":"123F8F39B626491983530EE5026003C2","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"可以看到，对于每一个参数采样，都生成了一个对应的后验预测值 Y。  \n\n我们可以进一步进行绘图："},{"cell_type":"code","metadata":{"collapsed":false,"id":"669553D25146458FA3F99F3D75D8B28C","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"sns.histplot(posterior_pred['Y'], stat='density', color='C0', label ='Posterior Predictive')\nplt.vlines(observed_data, 0, 0.16, color='r', label='Observed Data')\n\nplt.legend()\nplt.title('Posterior Predictive Distribution')\nplt.xlabel('Number of Students')\nplt.ylabel('Density')","outputs":[{"data":{"text/plain":["Text(0, 0.5, 'Density')"]},"execution_count":18,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/57AD8925A1E642F5A87686D8B4BC6620/s38bggowhx.png\">"],"text/plain":["<Figure size 640x480 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":18},{"cell_type":"markdown","metadata":{"id":"6212227DC427499FA70A11B047466F8A","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"## Summary and discussion  \n\n本节课我们学习了如何利用精确或近似的后验模型，对未知参数进行统计分析，包括：  \n1. 后验预测(均值和可信区间)  \n2. 基于后验的假设检验  \n3. 后验预测  \n\n😎到这里你或许发现了贝叶斯的独特支出，相比于之前学习的基于频率学派(frequentists)的统计方法。"},{"cell_type":"markdown","metadata":{"id":"35DD3915771542B0B2F013F3B2BE200A","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"贝叶斯分析结果的易解释性 -- 样本和总体的关系  \n\n- 贝叶斯分析根据观察到的数据 Y 来评估未知参数 π 的不确定性。  \n    \n    - 例如，在 100 名学生中，有 14 名学生在11点前睡觉  \n    - 根据这一观察结果，我们确定在11点前睡觉的学生比例 π 低于 0.20 的后验概率为 84.9%  \n        \n- 在频数分析中，往往根据对于总体性质π的假定值来评估观测样本Y的不确定性。  \n    \n    - 频率学派假设只要无限(long-run)的进行抽样，样本分布的参数就会近似总体分布。  \n    - 然后，根据零假设显著性检验(Null Hypothesis Significance Testing)和 p 值来判断样本和总体的关系。  \n    - 例如，频率学派假设 π 为 0.20时，我们观察到的样本(100 名学生)中，最多只有 Y = 14 (8%) 的学生在11点前睡觉。  \n        \n\n假设检验的哲学基础是波普尔的科学可证伪性  \n  - 它所回答的问题对人脑来说却**不那么自然**：因为我们实际上观察到了数据，但却不知道π，所以要解释一个假设相反的计算结果，可能会让人费解。  \n  - 然而，在检验假设时，更自然的做法是问 \"我的假设有多大可能性？(贝叶斯概率的答案），而不是 \"如果我的假设不成立，我的数据的可能性有多大？(频繁概率的答案）。  \n  - 鉴于P值经常被曲解，进而被误用，整个频繁主义和贝叶斯理论都越来越不强调P值。"},{"cell_type":"markdown","metadata":{"id":"8C6367AC3E4444BF994C939DCDA367EE","jupyter":{},"notebookId":"654741f6d74b63fed5dd17ca","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"两个学派的差异对比  \n\n频率学派  \n* 统计推断依赖于参数的**抽样分布**，即只要无限(long-run)的进行抽样，样本分布的参数就会有某种分布形式；  \n* NHST中的$p$值和置信区间的解读均依赖于“无限次抽样”的预设；  \n* 实际操作中，我们往往只会收集一次数据，并不会反复的进行抽样；有些情境中，预设“无限次重复抽样并不合理；  \n\n贝叶斯学派  \n* 假定参数本身是分布，不确定性一起存在于推断之中；  \n* 直接根据数据对先验信念进行更新；  \n\n|                     | 频率学派   | 贝叶斯学派   |  \n| ------------------- | ---------- | ------------ |  \n| 世界真相 (参数) | 固定       | 变化         |  \n| 概率                | 抽样的噪音 | 信念         |  \n| 推断过程            | NHST       | 贝叶斯定理   |  \n| 数据                | 存在噪音   | 固定         |  \n| 推断可更新性        | 否         | 是           |  \n| 主观性              | 前提预设   | 通过先验设定 |  \n"}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python","nbconvert_exporter":"python","file_extension":".py","version":"3.5.2","pygments_lexer":"ipython3"}},"nbformat":4,"nbformat_minor":2}